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DEVELOPMENT AND EVALUATION OF A KALMAN-FILTER ALGORITHM FOR 
TERMINAL-AREA NAVIGATION USING SENSORS OF MODERATE ACCURACY 
Gerd Kanning, Luigi S. Cicolani, and Stanley F. Schmidt* 
Ames Research Center 


SUMMARY 


The optimization and accuracy of a Kalman filter estimation algorithm for an 
integrated terminal area navigation system for passenger operations, using sensors 
and components representative of those expected to be commonly available on aircraft 
with instrument-flight-rules (IFR) and area navigation capabilities, are discussed. 
These sensors are body-mounted accelerometers and vertical and directional attitude 
gyroscopes, along with tactical air navigation aid (TACAN) or very high frequency 
omnidirectional radio range and distance measuring equipment (VOR/DME) , modular 
instrument landing system (MODILS) or microwave landing system (MLS), barometric and 
radar altimeters, and an airspeed sensor. This sensor set, available on several 
V/STOL aircraft at Ames Research Center is part of a digital flight control system 
called STOLAND. The principal investigative tool here is a simulation of the system, 
including sensors and their error processes; the discrete-time filter algorithm, 
together with the truncation errors and computational lags of the actual flight con- 
trol system; and a simplified model of the combined aircraft and control laws 
together with a reference trajectory command generator. Accuracy results are given 
as rms errors obtained from ensembles of 10 sample approaches along a reference tra- 
jectory. The filter optimization seeks to minimize computation time with negligible 
loss of accuracy and considers the appropriate selection of states, partitioning of 
the states into independent lower-order systems, and the minimum rate for processing 
navaid measurements to aid the acceleration measurements. Accuracy is investigated 
for the terminal area and for all filter states, including the basic (input) accuracy 
of measuring position and acceleration; the variation of estimation (output) accuracy 
throughout the terminal area, with maneuvering, location, flight direction, and axis, 
and its sensitivity to measurement accuracy; the trajectory dispersions and control 
activity excited by navigation errors; and a comparison of accuracy with that 
required to meet various terminal area safety criteria. It is found that estimation 
errors for this sensor set are nonstationary and nonisotropic in the terminal area; 
accuracy varies with maneuvering, and by an order of magnitude for the horizontal 
plane translational states during an approach, and differs by an order of magnitude 
with direction or axis at many points in the approach. Accuracy is sufficient for 
safety in IFR conventional and short takeoff and landing operations based on the use 
of VOR/DME and MLS. In addition, in automatic reference trajectory tracking, signif- 
icant tracking errors and control activity, excited by the estimation errors, affect 
ride quality and limit the usable control bandwidth, particularly for the horizontal 
plane motion. 


* Analytical Mechanics Associates, Inc., Mountain View, California. 



1 . INTRODUCTION 


This report considers the optimization and accuracy of a Kalman filter estima- 
tion algorithm for an integrated terminal area navigation system, using sensors and 
components representative of those expected to be commonly available on aircraft with 
instrument-flight-rules (IFR) and area navigation capabilities. The sensors assumed 
for this study are body-mounted accelerometers and vertical and directional attitude 
gyroscopes, which together form an inertial measurement unit (IMU) that measures 
inertial acceleration. A tactical air navigation aid (TACAN) or colocated TACAN and 
very-high-frequency omnidirectional radio range (VORTAC) , a modular instrument landing 
system (MODILS), barometric and radar altimeters, and an airspeed sensor are used to 
aid the IMU. 

This set of data types, together with a digital flight computer (Sperry 1819A) 
and other equipment, composes a digital flight control system (ref. 1) that is 
available in several vertical and short takeoff and landing (V/STOL) research aircraft 
operated by Ames Research Center (ARC) for flight test studies of estimation, guid- 
ance, and control algorithms. Although TACAN is a military system, it is functionally 
equivalent to the VOR and distance-measuring-equipment (DME) data of civil operations, 
and its bias statistics were set in this study to correspond to the less accurate 
civil data. MODILS is functionally equivalent to and represents the microwave landing 
system (MLS) facilities which are expected to support precision area navigation final 
approaches at many airports in the near future. 

These sensors, together with an appropriate algorithm, provide estimates of 
position, velocity, wind, and Euler attitude angles for use by the flight control 
logic. They differ from those used in earlier applications of Kalman filtering to 
area navigation (refs. 2-5) in that acceleration measurements are an order of magni- 
tude less accurate than the inertial grade IMU's of the earlier work, as a result of 
the fixed accelerometer errors and the maneuver-dependent dynamics of the pendulous 
attitude gyros. This system relies on intensive aiding of its IMU with position and 
airspeed data in order to obtain usable accuracy. 

The principal objectives of the study reported here were (1) to develop a 
Kalman filter algorithm for flight studies that yields the maximum estimation accu- 
racy inherent in the set of data types, within the constraints on computational 
requirements imposed by the airborne computer and flight control considerations; and 
(2) to evaluate the accuracy of the algorithm, particularly as it relates to the 
accuracy needed to support advanced, automatic, reference trajectory tracking 
operations . 

Estimation systems are needed to support IFR and area navigation operations in 
the terminal area. These operations are becoming increasingly common and are advo- 
cated at various levels of complexity (two-, three-, and four-dimensional area navi- 
gation). They are considered a requirement in the integration of V/STOL operations 
with existing conventional-takeof f-and-landing (CTOL) operations, and are generally 
advocated for all commercial operations by pilots, airline operators, and the govern- 
ment (refs. 6-8). This broad support derives from various proved or anticipated 
benefits, including increased safety, economy, and operational capability, and sig- 
nificant reductions in fuel, noise, and pilot and controller workload. 

The most demanding level of estimation and control system performance required 
in these operations is for automatic, four-dimensional (4-D) reference trajectory 
tracking and landing. The set of sensors studied here is expected to have marginal 
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performance for 4-D operations, but has the advantage of using equipment that is much 
less expensive than that used in state-of-the-art inertial grade IMU’s. Therefore, 
it is of interest to determine the maximum performance obtainable from these sensors 
and its suitability for 4-D operations. 

For this purpose, the optimal filter (minimum estimation error variances), when 
only accuracy is considered, is obtained from Kalman filter theory (ref. 9), provided 
system dynamics and measurement errors can be modeled with sufficient accuracy in the 
filter. However, the optimal filter is impractical to implement within the con- 
straints of computation storage and time limitations, because the number of variables 
whose uncertainties contribute to aircraft state estimation errors is unreasonably 
large. Thus, implementation of the Kalman filter always requires that only those 
variables be estimated that contribute significantly to performance in est ima ting the 
states required for the aircraft control. This entails deleting variables for which 
the information obtained in flight is negligible compared to the a priori information, 
including those variables that contribute negligibly to the estimation errors of the 
remaining states and those that contribute significantly to the errors but are none- 
theless poorly observable to the set of data types used. 

These and other approximations, such as linearization of the error state dynamics 
and the representation of measurement errors as Gaussian white noise, are ma de for 
tractability of the filter formulation, as well as to satisfy constraints on computa- 
tional requirements. Therefore, the objective of the filter development is to mini- 
mize computation requirements with negligible loss of estimation accuracy from the 
optimum. 

The principal investigative tool used in this study is a simulation of the sys- 
tem, including sensors and their error processes; the discrete-time filter algorithm, 
together with the truncation errors and computational lags of the actual flight 
control system; and a simplified model of the combined aircraft and control logic, 
together with a reference trajectory command generator. Accuracy results are given 
as rms errors obtained from ensembles of 10 sample approaches along a reference 
traj ectory . 

Various reasonable generic error models are used in the sensor simulations. For 
simplicity, the various types of off-nominal behavior observed in practice, other 
than data dropout, are neither included in the simulation nor are they within the 
scope of this study. These and other qualitative differences from the observed error 
processes exist, but the simulation suffices to determine filter performance trends 
with changes in the filter algorithm or in sensor accuracy and to evaluate the abso- 
lute performance of the system for nominal error behavior and statistics. 

The filter algorithm uses several implementation devices developed in the earlier 
work to solve various problems and to minimize computational requirements. These 
include (1) the square root formulation of the filter (refs. 10 and 11) to eliminate 
computational ill-conditioning and to ensure positive definiteness of the covariance; 
(2) measurement compression, to reduce the number and rate of scalar measurement pro- 
cessings by the filter, and (3) exponentially correlated random process models for 
states whose deterministic dynamics are unknown. In addition, this work attempts to 
settle several optimization issues aimed at minimizing the computation time required 
with only a negligible loss of accuracy. These issues include (1) the appropriate 
selection of state variables, (2) the minimum rate of executing the measurement 
processing computations, and (3) appropriate partitioning of the states into indepen- 
dent, lower-order systems. 


3 



A simplified generic model of the translational degrees of freedom of the combine! 
aircraft and control laws is derived. This model is independent of aircraft details 
and suffices to determine the trajectory tracking errors and control activity (mea- 
sured as the corrective accelerations required for trajectory regulation) excited by 
navigation errors. The model yields both analytical results for these relationships 
and simulation results for the effects of the present system's navigation errors on 
trajectory tracking performance. 

Sensor simulation models are described in the first section, and the filter 
algorithm is described in the second section. The remaining sections present a 
detailed analysis of the estimation accuracy achieved for all states throughout the 
terminal area and its sensitivity to measurement accuracy, the trajectory dispersions 
and control activity excited by the navigation errors, and a comparison of accuracy 
with that required for various terminal area operations. The available accuracy cri- 
teria generally reflect safety considerations within existing traffic separation 
standards and runway dimensions; they depend principally on the low frequency content 
of the navigation errors. Additional criteria for ride quality and control activity, 
which depend principally on navigation errors at or above the control bandwidth, 'are 
less developed than the safety criteria and were not considered here; however, they 
would significantly affect the suitability of an estimation system for use in auto- 
matic trajectory tracking. 

An analysis of the basic accuracy of the sensors in measuring runway referenced 
position is presented in appendix A, a model of the dynamics of the attitude gyro- 
scopes is provided in appendix B, and sensor accuracy in measuring runway referenced 
acceleration is discussed in appendix C. Appendix D contains an analysis of the 
effects of estimation errors on trajectory tracking errors. 


2. SENSOR MODELS AND MEASUREMENT ACCURACY 


Simulation models for all the sensors considered in this study are defined 
below; they are TACAN or VORTAC range and bearing; MODILS range, azimuth and eleva- 
tion; barometric and radar altimeters; vertical and directional attitude gyroscopes; 
and three-axis, body-mounted accelerometers. Generic models are given with parameter 
values, transmitter locations, and other details selected to correspond to the 
instrumentation and flight computer available at ARC and at its terminal area STOL 
test facility at the Navy Auxiliary Landing Field (NALF) , Crows Landing, Calif. 

TACAN and MODILS transmitting antennae are located, with respect to the runway, as 
shown in the view of the terminal area in figure 2.1. The locations of the antennae 
are appropriate for approach and landing operations and will permit realistic evalua- 
tion of the accuracy achievable with these sensors along STOL approach paths. The 
reference trajectory in this study is also shown in figure 2.1 and parameters defining 
this trajectory are noted in table 2.1 for later reference. The STOL runway is actu- 
ally painted on a longer standard runway so that the MODILS azimuth site is farther 
from the landing zone than is optimum for lateral position accuracy at landing; how- 
ever, this is readily accounted for in the evaluation. 


Reference Frames and Transformations 

The analysis makes use of several orthogonal coordinate frames; these are runway, 
path, body, and level heading axes. Runway axes, denoted (i , , k r ) (see sketch A), 

comprise a local vertical frame with the x-axis along the runway centerline; the 
navigation logic is formulated in this frame, and it is represented as an inertial 
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Figure 2.1.- Terminal-area navaid sites and reference approach path. 


TABLE 2.1.- TERMINAL AREA-NAVAID SITES AND REFERENCE APPROACH TRAJECTORY 


Leg 

Initial 

time, 

-sec 

Initial position, -km 

Initial velocity 

Acceleration 

X 

y 

z 

V, 

'knots 

K> 

-deg 

Y > 
-deg 

V, 

'g 

R c ’ 

-km 

1. Straight 

0 

0.629 

45.732 

-1.162 

140 

ESI 

0 

0 

-1.524 

2. Helix (it /2 ) 


0.629 

-.914 

-1.162 

140 

in 

-3 

0 


3. Straight 


-.895 

-2.438 

-1.036 

140 

HI 

0 

0.02 


4. Turn (ir/2) 


-4.298 

-2.438 

-1.036 

120 

180 

0 

0 

-1.219 

5. Turn (ir/2) 


-5.517 

-1.219 

-1.036 

120 

90 


-0.035 

-1.219 

6. Descent 


-4.298 

0 

-1.036 

96 

0 

-6 

-0.020 

OO 

7. Helix (2 tt ) 

828.7 

-2.807 

0 

-.883 

83 

0 

-7.5 

-0.010 


8. Glide slope 

929.6 

-2.807 

0 

-.381 

65 

0 

-7.5 

0 

OO 

9. Flare 

1014.9 

.024 

0 

-.008 

65 

0 

-1 

0 



5 










frame in the motion simulation and throughout this 
work. The atuomatic control is formulated in path 
axes for independent regulation of tracking errors 
along and normal to the path. These axes, denoted 
(ip, j_ p , kp) in the sketch, are oriented along the 
aircraft velocity vector (longitudinal axis) and 
normal to this in the horizontal plane (lateral, 
axis) and vertical plane (normal axis). Much of 
the estimation accuracy evaluation is carried out 
in this frame. Body axes, (ijj, jb> hb)» are 
required in the discussion of the IMU, and level 
heading axes, (i^, it. , k T ) , are also useful in 
some discussions under the small angle conditions 
of passenger operations. It is a local vertical 
system with i^ along the projection of ijj in 
the horizontal plane. 

Several transformations among these reference frames are used in the analysis; 
they are listed in table 2.2. The transformation of a vector's runway axis coordi- 
nates to its body axis coordinates is given by 

T br = E 1 (^)E 2 (0)E 3 (^) (2.1) 

where 4> , 0, and ^ are the usual Euler angles. The abbreviated notation, E-j^cf), is 
the transformation between orthogonal frames related by a single rotation taken about 
the ith axis and is defined in table 2.2. The transformation to path axes coordi— 
nates is 



T pr = e 2 ( Y )e 3 ( V (2-2) 

where y, ip v are the direction angles of the aircraft velocity vector (see sketch B) . 
The transformation to path axes based on the air velocity vector is obtained using 

the direction angles of the air velocity vector, 

T a » tv a - 

The transformation from air velocity path axes to 
body axes can be given from the above results as 

T bp - V^r ‘ M*)E 2 < e ) E 3 «Eh* Va >EhV (2 ' 3) 

a a 

Alternatively, this transformation can be given as 

T = E 2 (a)E 3 (-B)E 1 (* v ) (2.4) 

F a v 

where <J>y is the roll angle measured about the air velocity vector, and a and 6 are 
the usual angles of attack and sideslip, respectively, which locate V a relative to 
the body axes. Equations (2.3) and (2.4) can be equated to derive relations for any 
three angles in terms of the remaining five, as needed in the analysis. 

Last, the transformation from runway to level heading coordinates is simply 
E 3 0J») • 



VELOCITY VECTOR 
DIRECTION ANGLES 


Sketch B 


6 



TABLE 2.2- TRANSFORMATION MATRICES 


Single-axis rotations 


E,( a) = 


E 2 (o) 


E 3 (a) 


1 0 0 " 
0 cos a sin a 

0 -sin a cos cr_ 

cos a 0 -sin cf 

0 1 0 

sin a 0 cos 

“cos a sin a 0" 

-sin a cos a 0 

0 0 1 


Runway to body axes 






cos 0 

cos 

Y 







cos 0 

sin 

Y 




-sin 

0 


T u = 
br 

sin 

0 

sin 

0 

cos 

Y - 

cos 

4 > 

sin 

Y 

sin 

$ 

sin 

0 

sin 

Y + 

cos 

$ 

cos 

Y 

sin <I> 

cos 

0 


cos 

$ 

sin 

0 

cos 

Y + 

sin 

$ 

sin 

Y 

cos 

$ 

sin 

0 

sin 

Y - 

sin 

$ 

cos 

Y 

cos $ 

cos 

0 


Runway to path axes 


pr 


cos y cos Y cos y sin ¥ 


-sin Y . 


cos Y_. 


|_sin y cos Y sin y sin Y , 


-sin y 
0 

cos y 


Air-velocity path axes to body axes 


bp. 


cos a cos B -cos a sin B cos 4> v +sin a sin 4> v 


sin B 


cos B cos <}) 


|_sin a cos B -sin a sin B cos 4> v - cos a sin 4> v 
Runway to level-heading axes 

cos Y sin Y 


-cos a sin B sin d) - sin a cos d> 
v r v 

cos B sin <j> 


sin a sin B sin 4> v + cos a cos <f> v 


Lr 


i-sin Y cos Y 0 
L 0 0 1J 
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Navigation Aid Measurement Models 


Simulation models for the outputs of the sensors that aid the IMU are defined in 
table 2.3 and figure 2.2. For these sensors, a scalar measurement, Y, can be repre- 
sented in general as a deterministic function of the translational states and other 
variables, H(w), with measurement errors, Y, superposed: 

Y = H(w) + Y (2.5) 


Measurement errors arise from numerous independent sources in the equipment 
(sensors, receivers, transmitters, A/D devices) and from a priori calibration errors. 
Only the dominant error characteristics within the domain of use for a given sensor 
need be simulated for that sensor. Errors are represented as a sum of independent 
random processes with distinct statistical properties. 


where 


and 1 


N 

^ - E t ± + t 

i=i 


( 2 . 6 ) 


i. = (v . - £ . ) /t . 

^X v X 1 1 

v. - N(0,a. ) 

1 x 


t = Y - T[Y/T] - sign(Y)T/2 


Here, are Gaussian error components with zero means, standard deviations { a-j_ } , 

correlation times {t-^} which usually differ by one or more orders of magnitude. If 
is much smaller than the sampling interval of the estimator (0.1 sec here) then 
will be essentially independent from sample to sample. If x ^ is much larger 
than the terminal area flight duration, then £i will appear to be a bias for the 
flight. In this study, the Gaussian errors are modeled simply as a sum of the bias 
and moderately correlated sample errors. Last, t is a truncation error resulting 
from the limited resolution, T, with which data are represented in the digital com- 
puter. This error is approximately a sawtooth function in time. The value of T in 
the present context is given by the scaling selected for each data type in the 18-bit, 
fixed-point STOLAND flight computer. Assuming that t is uniformly distributed on 
[-T/2,T/2], the variance of the combined measurement errors is 


o 


2 



2 ^ 

a l + 12 


(2.7) 


Equation (2.6) defines the nominal error processes. In addition, there are both 
predictable and unpredictable situations in which a measurement is unavailable or in 
which its accuracy is significantly degraded from the nominal accuracy. The predict- 
able events result from such factors as scan limits of the transmitting antenna, 
signal strength loss with distance and in the fringes of coverage, shadowing of the 
receiving antenna, and dynamic instrument lags. The estimation algorithm imposes 
data admissiblity conditions of the general form 


1 The bracket [( )] denotes the truncation of the magnitude of ( ) to the nearest 
lower integer. 
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TABLE 2.3- SIMULATION MODELS FOR NAVAID AND AIR-DATA MEASUREMENT S a 





Error-model parameters 


Type 

Symbol 

Measurement function 
H(w) 

t = 10 b sec, 
0 

t = 1 sec, 
0 

T 

Error sum, 
a 

Measurement reception 
conditions, g(w) 

TACAN (VORTAC) 








Range 

Y 

tr 

|R - R t l 

61 m 
(305 m) 

30.5 m 

46 m 

70 m 
(307 m) 

0.3 < | R - R | < 800 km 



, /y t - y\ 






Bearing 

Y tb 

tan (x t - xj + *R 

0.57° 

(2°) 

0.1° 

0.25° 

0.58° 

(2°) 

0.3 < | R — R | < 800 km 
EL < 60° 

MODILS 








Range 

Y mr 

| R - R | 
■— -m 1 

6.1 m 

12.2 m 

18.5 m 

15 

0.3 < |R - R 1 < 16 km 
| Az | < 21.2' sm 

Azimuth 

Y 

ma 

tan_1 (y - y m ) /d xym 

0.17° 

0.07° 

0.1° 

0. 19° 

EL < 20° 
m 

Elevation 

Y 

me 

tan -1 (-z 1 /r 1 ) 

0.057° 

0.07° 

0.1° 

0.095° 

0.3 < |R - rJ < 16 km 

i_ — e i 

1.95° < EL < 16.5° 

|Az e | < 25° 

|-H < 20° 

Altimeters 








Barometric 

Y hb 

h R " 2 

30.5 m 

1.5 m 

0.075 m 

30.5 m 


Radar 

Y hr 

-z 

0 

0.6 m 

0.075 m 

0.6 m 

|<j>| <30°, 1 0 1 < 20°. h r < 60 m 

Airspeed 

Y Va 

lx - w| 

0 

0.6 mps 

0.003 mps 

0.6 mps 



Auxiliary quantities z lf r x> EL t , EL m » ^ are defined in figure 2.2. 


v o 






TACAN (VORTAC) 


'(x, y, 2 ) 


(x 


X 

t<yt' z t>/ a 


V&777777777777777. 


VERTICAL PLANE CONTAINING 
A/C AND TACAN 



MEASUREMENT FUNCTIONS AND 
AUXILIARY QUANTITIES 

R r = (x, y, z) T 
R tr = < x t - yt> z t )T 

(Ax t , Ay t , A^) = (x - x t , y - y t , z - zj 
EL t = tan ' 1 (-Az t A/ Ax t 2 + Ay t 2 ) 
H tr = Xx t 2 + Ay t 2 + Az t 2 
H tb = tan -1 (-Ay t /-Ax t ) + t// R 


MODILS 


(x, Y, z) 




VERTICAL PLANE CONTAINING 
A/C AND ELEVATION ANTENNA 


R m r ^ x irv ym ,z m^ 

(A x m » ^Y m ' Az m^ = " x nrr Y ~ Ym' z ” z m^ 

EL m = tan_1 < Az m /Ax m> 

^xym - \/ Ax m + Az m 

H mr"V /ax rr. 2 + 4 Vm 2+iz m’' ! 

H ma -t=n -UpS.) l Hma l<90- 
vu xzm 7 

Rg r = (xg, y e , Zg) 

(Ax e , Ay e , AZg) = (x - x e , y - y e , z - z e ) T 

Az e = tan" 1 (Ay e /Ax e ) 

z*| = Az e cos 5° - Ax e sin 5° 

^ = n /ax c 2 cos 2 5° + Ay e 2 + Az e 2 sin 2 5° 
Hme^an" 1 (^) lH ma K90° 


Figure 2.2.- Simulation models: VORTAC and MODILS measurement functions. 
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m = 1 , 2 , . 


M} 


( 2 . 8 ) 


{g m (w) > 0 


• 9 


which are conservatively designed to admit data only under conditions for which the 
nominal error models are known to be valid, except for unpredictable signal anomalies. 
These conditions define the effective domain of validity, W*, for each of the present 
sensor error models as 


W m = {w; 8 m (w) k 0}, m = 1, 2, . . . , M} (2.9) 

Thus, it is unnecessary to model off-nominal error behavior outside W*. Within W* 
some measurements, such as TACAN and MODILS, are randomly subject to signal dropout 
and various other anomalies; signal dropout is simulated for the study of navigation 
accuracy during dead reckoning, but otherwise the study of these anomalies and their 
effects on performance is beyond the scope of the present work. 

The measurement functions, error distribution parameter values, and measurement 
reception conditions of equations (2.5), (2.6), and (2.8) for the sensors of this 
study are listed in table 2.3 and figure 2.2; they are reviewed briefly below. 

Further discussion of navigation devices of the types used here can be found in 
references 12 and 13 along with additional bibliography; a summary of manufacturers' 
specifications for the STOLAND sensors is given in reference 14 and additional 
descriptive material on these devices, signal processing computations, and error 
analyses are presented in references 15-19. 

TACAN and VORTAC both provide measurements of the magnitude and bearing from 
magnetic north of a vector from the aircraft to the transmitter station (fig. 2.2). 
Measurement errors (ref. 18) are dominated by biases and imply large position fix 
errors compared with those of MODILS. Range and bearing can be received well beyond 
the terminal area. Bearing accuracy also degrades significantly at high elevations 
above the station and when passing near the station, a result of dynamic receiver 
lags excited by high bearing rates; these conditions are avoided or limited by exclud- 
ing bearing measurements at elevations above 60° and at distances from the station 
below 0.3 km. In addition, bearing bias can be subject to significant spatial varia- 
tions, depending on the multipath characteristics of the station environs. TACAN 
bearing is a military navaid which is functionally equivalent to the less accurate VOR 
navaid available for civil flight operations. Both types of navaids are used in net- 
works across the country in association with the National Airspace System. 

Error model parameter values are listed in table 2.3 for both TACAN and 
VORTAC instrumentation; these differ solely in bias magnitudes. The simulation 
results will be based on the civil-use VORTAC navaid. Accuracies for these 
navaids vary widely among stations and receivers (ref. 18) but only a single grade of 
equipment is represented here. The TACAN model corresponds to the facility in use at 
Ames Research Center; it represents an average-to-good station-receiver combination 
and is convenient in connection with local flight studies. The VORTAC model corre- 
sponds to low accuracy VOR and DME equipment; it will provide a somewhat conservative 
distribution of a priori biases for evaluating the accuracy encountered by a flight 
control system, but a distribution for which the effects of these biases and the 
possibilities for their in-flight calibration are more obvious in the Monte Carlo 
simulation testing. In this regard, the steady state estimation accuracy will be 
independent of the a priori biases to the extent that in-flight calibration is 
possible. 
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Some flight recordings and analysis of the TACAN station and equipment available 
at Ames is reported in reference 19. The observed errors are in reasonable agreement 
with the present model regarding noise, biases, and behavior near the station. 

Signal anomalies in the samples of reference 19 and as observed in other local flight 
experience include isolated large error events and lengthy periods of degraded high- 
frequency noise standard deviation, signal hangups, and dropout. Such events are 
largely outside the scope of this study, but are part of the practical estimation 
problem. 

MODILS is an experimental guidance system which provides range, azimuth, and 
elevation data (ref. 16). Antenna sites, coverage boundaries, and measurement geom- 
etry are indicated in figures 2.1 and 2.2 and table 2.3. Range and azimuth transmit- 
ting antennae are colocated nearly in the runway centerplane beyond the end of the 
runway; they provide measurements of the magnitude and angle from the centerplane of 
a vector from the antenna site to the aircraft (fig. 2.2). The elevation transmit- 
ting antenna is located a short distance laterally from the nominal touchdown point. 

It is a conical scanning antenna, tilted 5° above the ground plane and provides 
measurements of aircraft elevation above the ground plane. Azimuth/range and eleva- 
tion coverages are limited to volumes of the terminal airspace bounded in azimuth, 
elevation, and distance from the antennae as noted in table 2.3. These volumes 
suffice to cover the final portions of an approach from negative values of x. Ele- 
vation is further restricted to a narrow band of aircraft heading angles to exclude 
receiver shadowing effects during turns. 

Flight data (ref. 19) show gross agreement with the standard deviations of the 
simulation model. Sample signal error histories indicate some spatial dependence of 
biases, and the observed anomalies include infrequent and isolated, large elevation 
error events; there are also frequent episodes of linear divergence in range error 
(>250 m) caused by loss of receiver-transmitter synchronization. This latter type of 
anomaly is difficult to detect in the estimation logic. However, the poor reliability 
of the MODILS facility is not expected to be characteristic of future MLS’s and is 
only of special interest here. This is confirmed by flight experience with the 
Phase III Basic Narrow MLS now in use at the STOL test facility (ref. 20). The 
present MODILS simulation model is, therefore, only representative of such systems, 
not a fully realistic model of the observed MODILS errors. 

Altitude measurements are provided by the barometric and radar altimeters. 
Barometric altitude is obtained by sensing ambient free-stream pressure p a , and con- 
verting it to altitude from a stored, standard pressure-altitude model of the atmo- 
sphere, h (p a ) (ref. 15). On any given flight the model is calibrated to the mea- 
sured pressure at the altitude of the destination runway (Pr^r) > which is broadcast 
to arriving aircraft, so that barometric altitude is obtained as 

\ - h R + t h *<p a > - h *<p R >] 

The observed flight errors in this measurement are dominated by biases that vary with 
altitude above the point of calibration (because of departures of the actual pressure 
and temperature relations with altitude from the standard ones), and with speed and 
attitude (because of sensor calibration errors). A simplified model with fixed, 30-m 
rms bias is adopted for this study, but the effects of actual bias variations are 
noted in the evaluation where significant. Additional error types that need not be 
modeled in the present context are dynamic sensor lags and the large errors that can 
occur during flare and landing as a result of ground effects, as well as acceleration 
dependent errors excited by pitch up. 
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The radar altimeter measures altitude above the terrain based on the time 
required for the return of the transmitted radio signals. For commercial operations, 
it is used principally for final approach and landing. Here, the use of the radio 
altimeter is restricted to barometric altitudes under 60 m above the runway, where it 
can be assumed in this context that the terrain height is known. Typically, radio 
altimeters have very little bias error, but the random noise component increases with 
altitude, terrain roughness, and aircraft attitude angles; nevertheless, the measure- 
ments are much more accurate than those obtained with barometric altimeters for the 
restricted, low altitude range of its present use. A representative error model with 
null bias and fixed noise-standard-deviation of 0.6 m is adopted. 

The airspeed measurement (ref. 15) is derived from measurements of differential 
pressure, static pressure, and stagnation temperature, using standard atmosphere tem- 
perature and density models. A simple error model with no bias and a fixed noise- 
standard-deviation of 0.6 m/sec is used. 


Position Fix Accuracy 

Aircraft position can be calculated using three simultaneous measurements of 
independent functions of position; for example, TACAN range and bearing, and baro- 
metric altitude or MODILS range, azimuth, and elevation. This is done, for example, 
in the complementary filter navigation algorithm given in reference 14 to convert the 
actual measurements into equivalent measurements of the aircraft coordinates before 
processing the data. Position fix accuracy is also of interest here as the measure- 
ment accuracy against which any position accuracy improvements obtained from the 
filter can be compared. This topic is discussed in appendix A; there, general for- 
mulas for position fix accuracy are derived and applied to the combinations of posi- 
tion navaids of interest in this study to map the accuracy available in the ter mi nal 
area. 


Accelerometer Measurements 


A generic measurement model for the three— axis body— mounted accelerometers is 
given in table 2.4; the data, fm^,, measures the body axis components of specific 
force, f^, corrupted by four types of accelerometer output errors and a truncation 
error, t^, resulting from the selected resolution limits T^, with which these data 
are represented in the flight computer. The output errors are due to scale factor 
error, output bias and noise, and to misalignment of the accelerometer axes and body 
axes. The misalignment errors can be formulated as follows: the accelerometer axes 

{uj[} are located with respect to body axes by the angles (a i ,a i ) (see sketch C) where 
{a^} are the cone angles, , i^, Zu 2 , jj-, > Zu 3 , k^}, and is the clock angle 

which locates Uj[ on the cone. Nominal values of 0.25° are assigned to the cone 
angles and {o^} are uniformly distributed on (-180°, 180°). The measured body axes 
components are therefore related to the actual body axes components by the misalign- 
ment matrix, M, defined in table 2.4 (for small angles, {a*}). 


An expression for the resulting measurement errors is readily derived from the 
measurement model: 

s. 


fm b = f b ' fm b = 


a 2 sin a 2 


i_a 3 cos a 3 


a 1 cos a x 


a 3 sin a 3 


a x sxn o. l 


a 2 cos a 2 


*3 J 




Ub. 


5 b b b C t 


( 2 . 10 ) 
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TABLE 2.4.- ACCELEROMETER PACKAGE MEASUREMENT MODEL 


Measured specific force 


where 


h> = 

\ 

1 

+ \] 

< Mf b + «b 

+ b b > + 

c b 


b = 

a b 

" &b 






- 

1 

a 1 cos 

a i a i 

sin 

a r 

M = 

a 2 

sin a 2 

1 

a 2 

cos 

a 2 


_ a 3 

cos a 3 

a 3 sin 

a 3 

1 



Error type 

Distribution 

Parameter values 

Scale factor, s^ 


Si = 0.01, i = 1, 2, 3 

Axis misalignment 
Cone angle, 
Clock angle, a-^ 

~ N(0,a) 
a-L ~ U(-a,a) 

o = 0.25° 

a = 180° 

Bias, 

b ± ~ N(0 ,o) 

a = 0.015 g 

Noise: = (Vi - E,±) h 

Vi ~ N(0,o) 

(x,o) = (1 sec, 0.003 g) 

Digital resolution, T\ 


{T i } = {0.001 g, 0.001 g, i 








The corresponding error variances are 



Noting the parameter values listed in table 2.4 and the general restrictions on 
maneuver accelerations in passenger operations. 




= 1 g 


it can be calculated that the accelerometer biases dominate the measurement errors 
with some normal axis contribution from scale factor errors; that is, 

o 1 = 0.015 g 

a 2 = 0.015 g 

a 3 = 0.018 g 

The accelerometer errors are thus modeled principally as biases or as slowly varying 
errors of the order of 0.015 g to 0.02 g, with noise nearly an order of magnitude 
smaller. 


Measurement Error Models of the Vertical and Directional Gyroscopes 

The vertical gyroscope measures pitch and roll angle as the gimbal angles of a 
pendulous two-degree-of-f reedom gyroscope with its fixed axis along the body longi- 
tudinal direction and its spin axis controlled to track the apparent local vertical. 
The directional gyroscope measures heading as a gimbal angle of a two-degree-of- 
freedom gyroscope with its fixed axis along the body normal axis and its spin axis 
controlled to track the apparent magnetic north (ref. 12). The measured angles are 
represented as 


*g = * + *8 
6 8 - 8 + 5 8 
*s - * + h 
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The dominant errors are dynamic tracking errors in response to aircraft maneuvering; 
these are deterministic functions of aircraft acceleration and attitude histories and 
are generated from simulation models of the dynamics of the forced gyroscope. The 
simulation models for STOLAND system's attitude gyroscopes were obtained from con- 
tractor's notes and are given in detail in appendix B; response to control torques, 
fixed and friction drifts, and rotation of the local vertical with aircraft motion 
and Earth’s daily rotation are included. 

The mathematical description of these errors is complex but their behavior in 
the flight conditions and maneuvers of interest is readily described. First, in an 
extended, unperturbed, static equilibrium flight segment, errors relax to some steady 
value within ±0.25°. This steady state accuracy is the best-case gyroscope perfor- 
mance and is limited by background torques from friction and inertial rotation of the 
local vertical reference frame. The maximum gimbal rates used to drive the spin axes 
to their reference directions are as follows: 

Maximum gimbal rates 

Vertical gyroscope, both gimbals 0.033°/sec 

Directional gyroscope, leveling gimbal 0.050°/sec (2.11) 

Directional gyroscope, heading gimbal 0.110°/sec 

These values are an order of magnitude larger than needed to balance background 
torques, but they also limit the rates of reducing errors induced by prior maneuver- 
ing, as is illustrated by the transient response to initial errors shown in 
figure 2 . 3(a) . 

Second, longitudinal accelerations affect principally the pitch measurement 
error. These accelerations are small and within ±0.15 g for passenger operations; 
they can be constant during changes of aircraft reference speed or can vary stochas- 
tically as a result of control activity to track the reference trajectory. The ver- 
tical gyroscope acts to align the spin axis with the apparent vertical, _g - a., by 
nulling the sensed accelerations perpendicular to the spin axis. However, the control 
is cut off if the sensed acceleration magnitude exceeds a design value; this value is 
0.05 g for the pitch gimbal and corresponds to an angle of 3° between the spin axis 
and apparent vertical. The cutoff is in the middle of the normal range of longitudi- 
nal acceleration activity so that various distinct types of pitch error histories can 
result. A step change of longitudinal acceleration from steady state static equilib- 
rium conditions may be below or above the cutoff. If below, the pitch error rises 
nearly linearly to the angle between the true and apparent verticals at about the 
rate noted in equation (2.11). Acceleration magnitudes above the cutoff are usually 
developed rapidly by the aircraft before significant spin-axis misalignment from the 
local vertical can occur owing to control torques; after cutoff, the gyroscope is 
free in pitch (up to a time limit of 3 min) and subsequent misalignment develops very 
slowly, because of background torques. Both of these types of error histories are 
illustrated in figure 2.3(b). For acceleration magnitudes just above the cutoff, 0g 
reaches -0.5° in 100 sec and subsequently continues a slow drift. For acceleration 
magnitudes just below the cutoff, 0 is driven close to its steady state hangoff at 
-2.7° in 100 sec. For errors of this size, the pitch gimbal control would become 
locked out at a subsequent return to static equilibrium but this is countered by the 
cutoff time limit. 
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Stochastic longitudinal control activ- 
ity can delay the pitch error settling 
transient shown above; if activity is such 
that the sensed acceleration exceeds the 
cutoff for a significant percent of the 
time, then the settling transient is 
delayed and errors show some increase in 
steady state* An illustration of this 
effect was obtained by expanding the air- 
craft simulation to include a simple model 
of the control (described later); it is 
included in figure 2.3(b). 

Third, pitchover and flare are short 
duration maneuvers (not shown) , with very 
little acceleration orthogonal to the local 
vertical; therefore, they have almost no 
effect on 4>g, i/ig and only a minor tran- 
sient effect on 0g. 

Fourth, steady turns result princi- 
pally in large heading measurement errors, 
which are approximately sinusoidal and can 
rise to 4° or more (fig. 2.3(c)). The direc- 
tional gyroscope is controlled to track a 
sensed magnetic north; during turning 
flight, a large deterministic sinusoidal 
error in sensing north occurs and results 
in a corresponding control torque history 
which is the principal source of the head- 
ing error history seen in figure 2.3(c). 

For the vertical gyroscope, the turn excites 
no substantial errors if the gyroscope is 
initially well aligned; this is illustrated 
by the small sinusoidal errors, 0g, <J)g, of 
the order of 0.25°, seen in figure 2.3(c). 
The principal misalignment between true and 
apparent vertical in the turn is in roll 
angle, but the gyroscope’s roll gimbal con- 
trol is cut off if the sensed lateral 



(a) BEHAVIOR IN STATIC EQUILIBRIUM 



(b) PITCH ERROR FOR STEADY LONGITUDINAL 
DECELERATION 



(c) BEHAVIOR IN A STEADY TURN 

Figure 2.3.- Attitude-measurement errors. 


acceleration exceeds 0.1 g (this corresponds to a roll misalignment of 6° between 
apparent vertical and the gyroscope's spin axis). For most turns, lateral accelera- 
tion is well in excess of the cutoff and is developed rapidly during turn entry, 
after which the vertical gyroscope is free about the roll gimbal axis, and roll gimbal 
errors can change only slowly with the background torques. As in the case of longi- 
tudinal accelerations, it is possible to turn at accelerations just below the roll 
gimbal cutoff and eventually reach errors close to 6°. 


The altitude measurement errors for the last 500 sec of the test approach path of 
this study are shown in figure 2,4 along with time histories of the path axis compo- 
nents of trajectory accelerations. The attitude error behavior is consistent with 
that discussed above for isolated maneuvers; pitch errors are within ±1° and are 
associated principally with periods of acceleration below the cutoff value, 0.05 g 
(e.g., leg 3), and show sinusoidal behavior during turns. Roll errors rise to the 
order of 1° during the second and subsequent turns (legs 4, 5, and 7) as a result of 
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the initial misalignment of the ver- 
tical gyroscope at the turn entry 
(t = 730 sec). Heading error is 
within ±5°, and significant activity 
is confined principally to turns. 

All three errors show roughly con- 
tinual transient activity because 
of aircraft maneuvering during this 
part of the approach, and it is only 
on the glide slope (t > 925 sec) 
that the errors settle toward their 
steady state values for static 
equilibrium flight. Flight results 
with STOLAND gyroscopes (ref. 21) 
show qualititative agreement with 
the present error model but show 
that even larger error extremes 
than those in figure 2.4, of the 
order of 2° to 3° in pitch and 
roll, can result from maneuvering. 


Errors in Measuring Runway Axes 
Components of Acceleration 

The estimator uses measure- 
ments of the runway axes components 
of acceleration, which are derived 
from the outputs of the accelerom- 
eter and attitude gyroscope, using 

am r = T rb ( VW fn b + g r (2 ’ 12) 

errors. 
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Figure 2.4.- Attitude angles and derived 
acceleration: sample-case measurement 


Errors in the derived measurements are related to gyroscope and accelerometer measure- 
ment errors by (see appendix C): 


am r = 


"sin 0 f + cos 0 sin \p f 

r 2 r 3 

-sin 0 f - cos 0 cos ^ f 

r, T r 0 


cos ^ f r 
sin ip f 


-f 


Icos 0(cos ip f - sin \p f ) -cos \p f„ - sin ip f 

r i r 2 r i r 2 


<f> 


V 


+ (4> , 0 ,^) fmh 


(2.13) 


The notations am and fm indicate the measured acceleration and specific force, and 
{f r ^} are the runway axes specific force components. 

The effects of errors in the attitude gyroscope are given by the first term in 
equation (2.13). It is of interest to examine the magnitude and orientation of these 
effects relative to the path. For this purpose, expressions for the level heading 
coordinates of the gyroscopic term are derived in appendix C: 
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-80 g - a 2 * g 

g$ g + a J g I + H0T ' S 


,a 2 <j>g - ai 6 g/ 


(2.14) 


where (a x , a 2 , a 3 ) are the level heading components of acceleration (longitudinal, 
lateral, and vertical) and have maximum magnitudes of about 0.15, 0.35, and 0.1 g, 
respectively, in passenger operations. The result assumes that 0 is a small angle 
and that |a 1 | # |a 2 | « g. As seen in equation (2.14), roll error results in lateral 
acceleration error (to the order of 0.015 g) and some vertical acceleration error 
during turns (to the order of 0.005 g) ; pitch error results in longitudinal accelera- 
tion errors (to the order of 0.015 g) and in smaller vertical acceleration errors 
during speed changes. Heading errors cause no first-order acceleration error in 
static equilibrium and otherwise result principally in longitudinal errors (to the 
order of 0.03 g) during turns. 

The gyroscopic error effects are also rearranged according to the direction of 
the acceleration error in equation (2.14). Longitudinal acceleration errors depend 
on 0 g and on during turns, and lateral errors depend on $ g principally, and, to 
first-order, vertical errors are independent of attitude errors except for small 
effects, well under 0.01 g, during turns. 

The effect of the accelerometer package error, fm^, is given by the second term 
in equation (2.13). This error is principally a bias vector with components of the 
order of 0.015 g (rms), as noted previously, and is nearly a constant vector when 
viewed in body axes. When mapped to runway axes, this error is also constant during 
straight line flight with fixed attitude; however, during turns, only the vertical 
component is approximately constant and the horizontal plane components vary 
sinusoidally with heading. 


In summary, accelerometer biases and deterministic maneuver-induced gyroscopic 
errors dominate the horizontal plane components of am r , and fixed accelerometer 
errors dominate the vertical component. The principal errors are therefore fixed or 
of low frequency compared with the measurement sampling rate. Two sample histories 
of the path axes components of am for the test approach path are included in 
figure 2.4. In one of these, the accelerometer errors were nulled and the resulting 
errors are due to errors in the attitude gyroscope; as expected, longitudinal error 
is proportional to 0 g and to ij) g during turns, the lateral error is proportional 
to $g, and the vertical error is nearly independent of attitude errors. The second 
sample history contains randomly sampled accelerometer biases, including a large one, 
of the order of 0.05 g in the vertical axis; the effects of these biases are additive 
with the attitude error effects and are similar in magnitude. Although these errors 
are larger than the desired accuracy for navigation and control, they can be esti- 
mated and compensated to some degree in flight by the estimation algorithm to be 
described in the next section. 


3. A KALMAN FILTER TRANSLATIONAL- STATE ESTIMATOR 


The object of this section is to define a Kalman filter, terminal area, transla- 
tion state estimation algorithm for a digital flight control system, using the set of 
data types previously described. Both output accuracy and computational efficiency 
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are of acute importance in view of (1) the marginal accuracy of the IMU and VORTAC 
compared with that thought necessary for good IFR and automatic control, and (2) the 
large execution time of the Kalman filter computations in a multifunction flight com- 
puter compared to the fraction of each 0.05-sec computer cycle available for estima- 
tion computations. Therefore, the design objective is to minimize computation time 
while maintaining performance near the maximum that can be realized from the given 
sensors. 


Overview 


Estimates of a number of variables and parameters are required by the aircraft 
control; of these, the ones considered in this study are 

{R r > V r , W r , a r , <f>, 0, *} (3.1) 


The basic approach to their estimation uses measurements of attitude angles and body 
axes acceleration components to determine the runway axes acceleration; the trajec- 
tory states are then obtained by integrating the equations of motion: 2 


am r = T rb (*.e,V>)fin b + g r ] 


V = am 
r r 


R = V 
r r 


(3.2) 


The wind is estimated by combining airspeed and attitude measurements with the veloc- 
ity estimate. Since errors are present in these measurements, as outlined in the 
previous section, estimation errors using equation (3.2) will grow with time. How- 
ever, these errors can be estimated and corrected, using the other data types which 
are available and functionally related to aircraft position and velocity. Kalman 
filter theory (ref. 9) provides a formal basis for deriving an algorithm for this 
purpose. Measurements are made at discrete times, and corrections to the estimated 
variables are computed in proportion to the difference between the measurements and 
their values predicted from the current state estimate (measurement residuals); this 
provides an optimum (minimum variance) correction of the estimate in accordance with 
the relative accuracy of the measured and predicted values. Further, since the filter 
works on errors from the current estimate, the linear filtering theory can be applied 
using equations linearized about the current estimated state. 

Special devices used in the application include the square root formulation of 
the filter equations, exponentially correlated stochastic process models for some 
error states, and data compression. The square root formulation provides increased 
precision in finite word length computations and reduces computational errors owing 
to ill-conditioning to insignificant levels in the present application. It also 
enforces positive definiteness of the error covariance matrix. Error states whose 
deterministic dynamics are unknown (winds, measurement biases) are modeled as expo- 
nentially correlated random processes, with the results that the filter's covariance 
for these states degrades toward realistic values (to the a priori accuracy) during 


2 Runway axes are approximated as inertial in this study of terminal area naviga- 
tion since the errors involved are negligible compared with the measurements errors. 
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periods when no measurements are processed, and that realistic filter gains are com- 
puted when new measurements are subsequently processed. Data compression reduces the 
computation time required for measurement processing by summing multiple measurements 
of the same data type taken over a short time interval into a single nearly equiva- 
lent scalar measurement which can be processed in place of the individual 
measurements . 


State Estimation Equations 


The state X is an n-vector of the variables to be estimated and is assumed to 
satisfy the differential equation 


X = f (t) 

from which the estimated state is computed as 


(3.3) 


I X(t k ) + JT X(T)dT t k < t < tktl 

k (3.4) 

X(t k+1 ) + t = t k+1 

^ ere * the times (t^, k=l, 2, . . .) are those discrete times at which a correction, 
AX, is provided by the Kalman filter algorithm; at these times, & is changed dis- 
cretely, and between these times X is obtained from measured or assumed values of X 
or its integral. Equations for the variables of interest are given in table 3.1. 

These include measurement biases _b, used in the filter algorithm in addition to the 
variables noted in equation (3.1). 


Estimation Error State 

Between filter corrections, the estimation error state X is assumed to be 
governed by linear perturbation equations forced by a Gaussian white noise vector 
process, r|(t): 


X(t) = X(t) - X(t) 

x = fx + n 

n ± ~ N (0 ,q ± ) i = 1, 2, . . . , n 

The general solution for X(t) is given from the transition matrix, <J>(t,t k ), which 
can be satisfactorily approximated by a truncated Taylor series for the short inter- 
vals that occur in this work: 


(3.5) 

(3.6) 


X(t) = <Kt,t k )x(t k ) + u(t) 


(3.7) 


where 

<Kt k + 6,t k ) = e" F(S = I - F6 + F 2 6 2 /2 + . . . 
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TABLE 3.1.- STATE ESTIMATE AND ERROR STATE VARIABLES 



Estimate between corrections, 
C k * * < t k+i 


Estimate at a correction 


X(t) - X(tt) + f XOOdT ^ (t k+i ) ^' l ' C k+l ^ + 4X(t k+l ) 


>t 6 r (t) * \ (t k> + /' V*> dT V<W + 4 VW 


'r (t > " V c k> + .f \ (1 


VW + 4 V 4 k«> 



a r (t) = am r (t) + Aam r (t k ) 


e s + 49 


W (t) = W (t. ) 
r 7 r k 


b_(t) = Mt^) 


am (t, , ) + Aam (t, , „ ) 
r v k+i r k+i 7 


+ A0 


\s+4 


w r (t k ) + ™ r (t k+1 ) 


b(t k ) + ib(t k+1 ) 


Error state variables 


position estimation error 


velocity estimation error 


acceleration measurement bias, runway axes components 
acceleration measurement bias, level heading components 
vertical acceleration measurement bias 




attitude gyroscope biases 


W ,W 
x y 


wind error, runway x,y-axes components 


b ,b , ,b u , VORTAC range and bearing, and baroaltimeter calibration 
tr tb hb 

errors 





<Kt >t k )n(T)dT 


u(t k + 6) 



The selection of the filter's error state variables affects both estimation 
accuracy and the time required for executing the filter equations. For maximum 
accuracy, the state should include not only errors for the variables required by the 
control but also errors for all variables and parameters of the measurement functions 
for which values must be given in order to define the relation between measurements 
and the control variables. However, the filter execution time rises with the cube of 
the number of state variables to be estimated (ref. 5); as a result, it is essential 
to omit variables that do not significantly affect the accuracy of the control varia- 
bles or that if included do not substantially improve the estimates of the control 
variables. These are variables for which the amount of on-line information is 
small compared with the amount of a priori information (they are poorly observable 
to the measurements). The effect of their a priori estimation errors on the 
accuracy of the remaining variables can range from negligible to significant. 


The process of choosing appropriate error state variables is specific to each 
application and ultimately rests on simulation testing. The variables considered in 
the present study are listed in table 3.1. The horizontal plane components of the 
wind are included but the vertical component is omitted, because the on-line informa- 
tion on it is negligible and its error has little effect on estimation accuracy for 
the inertial vertical axis motion. Measurement biases are included for VORTAC and 
baroaltimeter , since their corresponding position errors are statistically large 
a priori compared with the desired position accuracy; moreover, they can be estimated 
ou— line by the more accurate M0D1LS and, to a limited degree, by other data types. 
MODILS biases are nontrivial, but they are omitted for lack of independent calibrating 
measurements of sufficient accuracy among the navaids . Transmitter location and runway 
parameter errors have negligible effects and need not be considered. In addition, 
this study considers three candidate formulations of the acceleration measurement 
errors in an effort to maximize the filter's capacity to detect and compensate the 
significant low-frequency errors encountered during maneuvering. The first formula- 

represents these errors simply as biases of the measured runway axis components 
of acceleration and results in a filter with 14 states; the second uses attitude and 
vertical acceleration measurement biases (15 states), with the object of detecting 
the gyroscopic error transients more accurately and possibly improving the attitude 
accuracy; and the third is formulated as biases of the measured level-heading-axes 
components of acceleration (14 states), with the object of avoiding the sinusoidal 
inertial acceleration measurement errors that occur during turns for the runway axes 
formulation owing to the accelerometer package bias vector. 

State equations for the variables of these three filter cases are given in 
table 3.2, from which the matrices F,<f> of equations (3.6) and (3.7) can be given by 
inspection. The wind and measurement error state variables are all modeled in an 
ad hoc manner by first-order differential equations forced by white noise; that is, 
they are represented as exponentially correlated random processes. The variance and 
time-constant (a , t) for each variable are selected to correspond to the time scale 
at which significant changes occur and the steady-state accuracy reached in the 
absence of new navaid measurements. The forced term u is a random variable whose 
variance a [1 — exp(-2S/T)] is such that the desired steady-state variance is 
obtained (E[z 2 ] -> a 2 ). 
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EQUATIONS 


Parameter 

values 


sec 


0 


0 


-03g, .03g, . 03g 

; 

25, 25, 

1 

6.1 mps 

20 

6. 1 mps 

20 

05 m 

‘O' 

2° 

1 

! 

1 

i 


30.5 m 


10 ' 








Error State Covariance 


Measurement processing in the filter is referenced to discrete times (t k ) at 
which the error state covariance, P(t k ), defined by 


P(t k ) = E[X(t k )X T (t k )] (3.8) 

is required. Its propagation between these times follows from equations (3.7) 
and (3.8) as 


P(t k+i ) + Q 


(3.9) 


where 


Q = E[u(A)u T (A)] = diag[a?(l - 2A/x i )] 


In forming Q, noise components for the kinematic states and correlations between the 
components of u have been neglected. The resulting matrix is diagonal, with ele- 
ments and parameter values as given in table 3.2. Equation (3.9) describes the 
theoretical behavior of accuracy with time. A, in the absence of new navaid measure- 
ments. Both terms vary with A, and, in view of table 3.2, they combine such that 
the variance of each error state (diagonal terms of P(t k+1 )) either increases indef- 
initely with A (e.g., R r , V r ) or increases to the steady state value, o^, which 
reflects the^ appropriate sensor accuracy (e.g., am r , <ji) or the a priori accuracy 
(e.g., b tr , W x ). 

The square root covariance will be used in the filter formulation. This is an 
n x n matrix, W, such that P = WW . An n x 2n square root of the propagated 
covariance is readily given by separating equation (3.9) into the form 

P(t k+i } = ^ (3-10) 

where 


T 

and then the matrix A can be reduced to an n x n upper triangular, square root 
matrix, W(t k+1 ), using Householder’s algorithm, # (ref. 10), since the product, AA T , 
is invariant for these operations: 


w T(t k+ i) =^(A T ) 


that is, equation (3.10), together with Householder’s algorithm, is used to advance 
the square-root covariance in time. 
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Measurement Models 


The filter computes corrections, AX, by processing samples of the available data 
types, which it assumes can be modeled by a deterministic function, H, with a super- 
posed Gaussian white noise process, Y: 


= H (X,p ) + Y m = 1, 2, . . . , M 
m m o m 


(3.11) 


Here, M is the number of data types and p refers to those parameters on which the 
measurement depends, in addition to X, but for which a priori values p Q are 
assigned and estimation errors p are neglected. The Gaussian white errors only 
approximately represent the errors between H m (X,p 0 ) and the actual measurements; 
these can include, for example, truncation errors, moderately correlated sample 
errors, and nontrivial biases. The noise variance q m is, therefore, adjusted 
empirically from simulation tests or recorded flight data to obtain the best fit. 

The measurement models and parameter values for this study are listed in 
table 3.3. A comparison with the simulation models (table 2.3) indicates the differ- 
erences in error model details. In addition, note that the MODILS elevation measure- 
ment has been modified to a derived altitude measurement and that the airspeed 


TABLE 3.3- MEASUREMENT MODELS 


Y = H (X , p) + Y 

Y - N (0 ,q) 


Measurement 

Symbol 

Measurement function, H(X,p) 


TACAN (VORTAC) 


[(x - x t ) 2 + (y - y t ) 2 + (z - Zj .) 2 ] 1 / 2 + b tr 


Range 

Y tr 

91.4 m 



Jy t - y\ 


Bearing 

Y tb 

t£ln \ x t " x ] + ^ + btb 

1° 

MODILS 




DME 

Y 

mr 

[(x - x m ) 2 + (y - y m ) 2 + (z - V > 2 ] 1/2 

45.7 m 

Azimuth 

Y 

ma 

tan -1 (y - y m )/Kx - x m ) 2 + (y - y m ) 2 + (z - z m > 2 ] 1/2 

0.3° 

Elevation 

Y' 

-z 

0.005 v 1 m 


me 


(see note) 

Baro altimeter 

Y hb 

_Z + h R + b hb 

1.5 m 

Radar altimeter 

Y hr 

-z 

1.8m 

Air velocity 

Y ia 

x - W 

X 

1.8 m/sec 


Y- 

ya 

y - w 

1 . 8 m/sec 



y 



Note: r x = [ (Ax e cos 5° + Az e sin 5°) 2 + Ay 2 ] 1 / 2 ; Ax 0 = (x - x e > , etc. 
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measurement has been combined with heading to yield derived measurements of the hori- 
zontal plane coordinates of the air velocity vector. The derived altitude measurement 
Y me is introduced to facilitate division of the state variables and measurements into 
separate filters for the horizontal plane and vertical axis motions. This measurement 
is derived from the -elevation measurement and current position estimate, using an 
expression derived from the simulation model. 


Y' = z + Ax tan 5° 
me e e 


r, tan Y /cos 5° 
me 


This is modeled in the filter as 


(3.12) 


Y' = -z 
me 


+ Y' 
me 


Using the simulation model, the error Y me can be related to the elevation measure- 
ment and position estimation errors. After neglecting second-order effects, this 

yields the standard deviation of Y' as 

me 


which increases with distance from the elevation antenna. 


The derived air velocity components are calculated from airspeed and heading 
measurements, using 


Its measurement model is 



= Y 


Va 


( cos \p 
sin \p 


g 

g 



(3.13) 


The measurement error can be derived by expressing the air velocity components in 
terms of airspeed and heading, using 


Va 


r 



Va 


P 


(3.14) 


where 


Va r " ( VV*a )T 
Va p - (V a ,0,0) T 

The transformation Tp r can be rewritten as Tp^T^j. and evaluated using table 2.2, 
after which the horizontal plane coordinates can be obtained as 
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X 


= Va 


'cos (ip + A^) 


,sin (\p + A\p), 


(3.15) 


where 


Axp = 3 cos <p - a sin p + HOT's 


This result assumes that a, 8, and 0 are small angles. The measurement errors for 
the derived air velocity components can now be given as: 




r x 


'cos Ip' 


= Y. 


Va \ . . 

sin \p J 


+ Va(i - A\p) 

o 


^-sin iJjN 


cos \p; 


+ HOT’s (3.16) 


where the error sources (Yy a , \p R9 and A\p) are the airspeed measurement error,' the 
directional gyroscope error, ana the unmodeled effects of (a, 3). The two error vec- 
tors in equation (3.16) are mutually orthogonal and in the horizontal plane along the 
direction of the aircraft heading and lateral to that heading. Thus, for the 
derived measurements, the error in measuring the longitudinal air velocity component 
is given by the airspeed sensor error, Yy a , as expected, and the error in the lateral 
direction is proportional to both the error in the directional gyroscope and to 
unmodeled angle effects; the error in the lateral direction can be of the same magni- 
tude as the air velocity component in this direction. The complex errors of the 
derived measurements in equation (3.16) are represented only approximately in the 
filter; measurement error correlation for the two components is neglected, and 
standard deviations of 1.8 m/sec were assumed for these errors. 


Measurement Preprocessing 

Navaid measurements are received and processed at discrete times as outlined in 
sketch D. Measurements are received and accumulated at a rapid rate (10 Hz) and 

these are processed by the filter to compute a 
correction to the state estimate, AX, at a 
slower rate (rates of 0. 1 to 2 Hz are studied). 
For convenience, the preprocessing logic refers 
all measurements taken during (t^tk-K) to a 
reference time, t^. Further, all admissible 
measurements of a single type are accumulated 
as a single nearly equivalent measurement for 
the interval. The state estimate correction is 
then computed by sequential processing of a 
single scalar measurement of each type. 

The processing uses measurement residuals; that is, the difference between the 
measurement and its predicted value from the current state estimate: 

y (t) E Y (t) - H (x(t),p) m = 1, 2, . . ., M (3.17) 

m m m 

Its relation with the estimation error at the reference time is obtained by lineariz- 
ing equation (3.11) about X(t) and using equation (3.7): 


MEASUREMENT RECEPTION TIMES 

/ 

l 1 x 2 ■••• 

l ■ I i I ■ I i 1^ 1 i 1 ( 


r 

l k 


l k+1 


FILTER UPDATE TIMES 

Sketch D 
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where 


y m (t) - Vt)*(t,t k )x(t k ) + \ 


m = 1, 2, . 


M 


(3.18) 


• • > 


h „(t) - [^H(X,p o )]. 

"(O-lm) 

The forced solution, u, in equation (3.7) can be neglected here, provided the filter 
update interval is small relative to the correlation times of the estimation error 
variables . 


The measurement gradients (h m ) for the present application are listed in 
table 3.4. In this table, the state variables are separated into two sets associated 
horizontal plane and vertical axis motions, respectively, to facilitate the 
analysis of separate filters for these two sets of variables in a later section. The 
filter algorithm is implemented with this same ordering of the variables to facilitate 
simulation study of the separate filters. 


Degraded navaid data can occur randomly; this refers to data which may have no 
relation to the state or contain errors generated by processes significantly larger 
than postulated in the filter design. In the present context, this can occur in iso- 
lated samples or for periods of any duration, and can occur as signal dropouts, signal 
hangups or divergence, or increased noise variance. We seek to exclude such data 
which, if processed, can result variously in large amplitude excitation of the fil- 
ter’s impulse response, or divergence of the estimation error for a period of time or 
excessively noisy estimates. Much of the degraded data is excluded by receiver 
validity checks and conservatively computed coverage boundaries for admitting data 
(given in table 2.3). However, this does not suffice to exclude all degraded data, 
so the filter rejects residuals which are large compared to its standard deviation 
computed from the filter’s covariance; that is, 


If: 



h Ph T + 
m m 



then: delete y 

m 


(3.19) 


Values of 2 to 4 have proved satisfactory for c in empirical tests with real data. 
This device succeeds in excluding dropouts of any duration and isolated large residual 
samples which would otherwise excite the filter's impulse response. Note that large 
residuals also result from sufficiently poor prior estimates as well as degraded data 
so that normal measurements would be locked out if errors in the predicted measure- 
ment exceeded the threshhold. For example, if VORTAC or baroaltimeter biases for- 
tuitously exceeded the threshhold then the swith to one or more MODILS measurement 
can become impossible or require reinitialization. However, this is a low probability 
event and was never observed with real data or in simulation tests with the normal 
system model. 

No attempt is made here to treat signal hangups or divergence and increased noise 
levels; schemes for signal hangups are unknown and various schemes to detect and treat 
changes in noise model parameters are available (e.g. , some schemes and applications 
are discussed in refs. 22 and 23) but are elaborate and probably ineffective in the 
present application. 
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TABLE 3.4.- MEASUREMENT GRADIENTS AND RESIDUALS: 14-STATE FILTER 


1 

rr 


cos ip^. cos 

EL t 

sin ip cos EL t 

0 

0 

0 

0 

1 

0 

0 

0 

-sin EL t 


0 

0 

0 















1 






X 



sin ip. 


-cos tp 









1 






y 



L 


t 

0 

0 

0 

0 

0 

1 

0 

0 

1 o 


0 

0 

0 



1 V 






'tb 


d 


d 











X 


xyt 


xyt 









1 




















1 






y 

y mr 


cos Az cos 

EL m 

cos Az 

0 

0 

0 

0 

0 

0 

0 

0 

I cos Az sin 
1 

EL m 

m 

0 

0 

0 


X 



cos Az cos 

EL 

cos EL 









sin Az sin 

EL 

m 





9 




m 

m 

0 

0 

n 

0 

0 

0 

0 

0 

I 

0 

0 

0 



^ma 


d 


d 

u 

1 d m 



b tr 


m 


m 









m 






~ 


= 












i 






b tb 

^xa 


0 


0 

1 

0 

0 

0 

0 

0 

-1 

0 

1 0 


0 

0 

0 


*x 














i 







y • 


0 


0 

0 

1 

0 

0 

0 

0 

0 

-1 

1 0 


0 

0 

0 


y 

ya 













1 







y 


0 


0 

0 

0 

0 

0 

0 

0 

0 

0 

! -i 


0 

0 

0 


z 

me 













i 




















i 






z 

y hb 


0 


0 

0 

0 

0 

0 

0 

0 

0 

0 

i-i 

i 


0 

0 

1 


z 

Jhr_ 


0 


0 

0 

0 

0 

0 

0 

0 

0 

0 

]-i 


0 

0 

0 


b hb 


3H 3H 3H 

Notes: In the 15-state filter -r-r = -rr = irr = 0 for all measurements. 

o<p du dip 


ip , EL , EL , d . d are defined in appendix A, table Al. 
r t t 9 m 5 xyt m 



Multiple measurements of a single data type received during the interval 
(tk’tk+i) can be compressed to a single scalar measurement that is nearly equivalent 
to the set of measurements by summing residuals, gradients, and variances: 


n. 


m 


y sm ^ y m^i^ 


1=1 

% 


sm 


1 sm 


1=1 


_ „ 1. 4 
q n 
^m m 


m = 1, 2, 


, M 


(3.20) 


where n m is the number of admitted measurements for the mth data type. In gen- 
eral, two sets of measurements, {(h^,qi), i = 1, 2, . . . , n} and 

{(hi^i)* i = 1» 2, . . . , n) are equivalent, that is, yield identical accuracy P 
after processing, provided they have the same information, that is, provided 


I 


T 

hTh. 

l l 


n 


I 


1=1 


h'. T h! 
i i 


(3.21) 


In the simple case in which all measurements are identical and denoted (h,q), the 
information can be written variously as 


n 


E 


1=1 


T T 

h7h. - T, h h 
ii h h s s 

= n = 

q n q 


where h s is the summed gradient, nh. Thus, { (h-j_ , qj_) } is equivalent to the single 
summed measurement (h s ,nq) which can be processed in place of { (hj^q^) } at a consid- 
erable savings in computation time and with no loss of accuracy. More generally, any 
collection of measurements {(h^q^)} can be compressed to a minimum equivalent set, 
{(hi>qi)} i n which (h^} is a basis of the space spanned by { h-^} (ref. 24). In the 
present application, the measurements of a single data type taken during (t k ,t k+1 ) 
are nearly identical and can be equivalenced to the single summed measurement of 
equation (3.20) with negligible error. In addition, the actual measurement errors 
can be moderately correlated at the 10-Hz sampling rate used here. Therefore, the 
summed measurement is weighted conservatively in equation (3.20) by increasing its 
error variance by the factor n 0,4 over its value in the case of independent errors, 
nq. This adjustment was selected empirically using recorded flight data. 


Measurement Processing Equations 

Each summed measurement is processed using the following square root, Kalman 
filter algorithm (Potter's algorithm, refs. 10 and 11). The variance of the residual 
s m , optimal (minimum variance) gain K, error state correction Ax, and posterior 
square root covariance W, are given by 
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T 

p = W h 


sm 


s = pp + q 
m M sm 


K = Mp/sm 

a = AX + K( y S m “ h sm AX > 

w T - w T - p k t /[ i + (q sm /s m ) l/2 ] 


>m = 1, 2, . . ., M (3.22) 


Here, the equal sign indicates replacement in the computational sense and all quanti- 
ties are associated with the reference time, t^. Note that in computing AX, the 
summed residual, y sm , modified to account for the error correction already accumu- 
lated but not used to calculate the residuals, y sm . 


Filter Initialization 


Starting values of the state estimate and estimation error covariance are 
required. In general, the initialization can be formulated as in equation (3.23). 

An initial set of measurements, or assigned a priori values, designated Y D , suffice 
to define uniquely the state variables; that is, and X are n-vectors and 
H(X,p 0 ) has a nonsingular Jacobian with respect to X. Then, the initial state 
estimate is given from the inverse function, H” 1 (Xo>Po) J an< * t * ie initial state esti- 
mation error, X Q , and its square-root covariance, W, are given from the Jacobian of 
the inverse function, J, and from the measurement errors, Y, and their variances, 
{q^}. The measurement errors can be assumed independent so that D is diagonal: 

Y - H(X,p ) + 1 


— o 


X = H -1 (Y ,p ) 
<"> — 0*0 


J = [V y H 1 ] y 


X = -JY 


o 
D = 

P = 
o 

w T = 


E[Y Y T ] 
-o-o 


J DJ 
•S J T 


diag{q^} 


(3.23) 


In the present application, the following simple initialization can be used at 
arrival in the terminal area. The winds, measurement biases, and vertical velocity 
are all assigned their a priori mean values (zero) and then the kinematic states are 
calculated from TACAN or VORTAC, baroaltimeter , airspeed, attitude, and accelerometer 
measurements. The corresponding initial estimate, expressions for the initial state 
estimation errors, and standard deviations for the initial measurement errors are 
given in table 3.5. Last, the rows of W£ are listed in table 3.6 for the 14-state 
filter; the ith row is associated with the ith measurement. 


32 



























and W Q can be formed with its rows taken from table 3.6 in any order. 

The initial estimate defined above can be in substantial error; for example, the 
initial measurement bias errors equal the a priori calibration errors, the initial 
horizontal plane velocity and wind errors both equal the actual wind, and the initial 
vertical velocity error equals the actual vertical velocity. After initialization, 
the filter converges to some M steady-state M covariance; this convergence is rapid and 
the covariance (theoretical accuracy) becomes independent of the initial covariance 
for those states for which new measurements provide information that rapidly out- 
weighs the initial information in Pq 1 . For other states, information is acquired 
very slowly or is unavailable until later in the approach and their accuracy is 
dominated by the a priori accuracy. 

In general, the estimation error can diverge following initialization for suffi- 
ciently large initial errors as a result of dynamic and measurement nonlinearities 
neglected in the filter. In some applications, convergence is sensitive to these 
nonlinearities and special devices are useful (e.g., ref. 25). In the present work 
the boundaries of convergence were not studied systematically, but we note that no 
divergence was ever encountered in simulation tests with 2a and 3a biases and winds 
and initialization at various distances from the runway. Thus, initial convergence 
of the estimate in the terminal area navigation problems appears insensitive to 
initial errors for the present initialization scheme and measurement error statistics. 


A Kalman Filter Algorithm 

A computational flow diagram and equation summary of the Kalman filter algorithm 
for this study are given in figure 3.1 and table 3.7. 

The computational flow is structured to permit reduction of the required compu- 
tation time in the flight computer with negligible loss of accuracy. For this pur- 
pose, computations are separated into parts corresponding to (1) integration of the 
equations of motion, (2) data reception and preprocessing, and (3) measurement pro- 
cessing with the Kalman filter to obtain a new estimate of the error state. This 
permits execution of these parts at different rates which can be separately optimized. 

Integration of the equations of motion requires little time and is executed at 
the maximum possible rate allowed by the computer’s cycle time (20 Hz) for maximum 
accuracy. Execution time for the measurement reception and preprocessing equations 
is moderate and is done at 10 Hz. More frequent samples have increasingly correlated 
errors and would not significantly increase the actual rate of accumulating informa- 
tion; less frequent sampling saves little computation time and loses significant 
information for some data types. In addition, the preprocessing deletes some data 
types whose contribution to estimation accuracy becomes negligible when more accurate 
and functionally equivalent data types are available; for example, VORTAC is deleted 
when MODILS is available. 

Execution time for the algorithm is dominated by the Kalman-f ilter measurement 
processing equations (ref. 5); therefore, the required computation time, as a percent 
of real time, depends principally on the rate of processing scalar measurements. 
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(1) NOTE: FOR 14 AND 15 STATE FORMULATIONS, 

THE FILTER OUTPUT IS RESTRICTED AS FOLLOWS: 

A<j> = A9 = A\p = O 14 STATE FILTER 
Aa r s(0, 0, Ai) 15STATE FILTER 


Figure 3. 1 Kalman-f ilter estimation algorithm: logical flow diagram 




TABLE 3.7.- KALMAN-FILTER ESTIMATION ALGORITHM: EQUATION SUMMARY 


Rate Equation Comments 

10 Hz Preprocessing logic: accumulate measurements for each data type, m = 1, 2, . . .,9 


If measurement is valid, continue 


Validity test 


y = Y - H (x) 

m m m 

If |y I < ca , continue 
1 m 1 m 


Residual 


Data rejection test 


^ s m ^ s m + 


»« - [ra m (i)]*(t,e k ) 

hs m h -m + 


Residual sum 


Measurement gradient with xCt^.) 


Gradient sum 


n = n + 1 
m m 


Number of accumulated 
measurements 


q s = q (n ) ] 
n m m 


Error variance of summed 
measurement 


1 Hz Process summed measurement for each data type selected and update estimated state 


If data type selected, continue 

T T ^ 

p = W h 

sm 
T ^ 

S = pp + q 
m n sm 

K = Wp/s 


Selection test 


Variance of summed residual 
Kalman gain 


AX = AX + K(y - h AX) 

; sm sm 7 

W T = W T - pK T /[l + (q /s ) l/2 ] 

sm m 

a = (s /n ) x / 2 
m m m 


m = 1 , . . 


Posterior error state estimate 

Posterior square-root covariance 

Standard deviation, single 
measurement residual 


(h sm’ ( Wy S in ) = (0 ’°’ 0) 
x = x + Kt k+i ,t k )Ax 


Reset measurement sums 

Posterior state estimate at 
current time, t, , , 


AX = 0 


w T = j^(a t ) 


Reset state correction 

Propagate square-root covariance 
to c k+i 

Upper triangular nxn/p, 
(Householder algorithm) 
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Maximum accuracy is achieved by processing all available data as soon as received, 
but this can be traded for computation time by reducing the processing rate. Factors 
in optimizing this rate are (1) the IMU accuracy, which governs the rate of accuracy 
loss for the kinematic states between filter updates; and (2) the measurement summing 
logic, which compresses measurements taken during the filter update interval into one 
measurement per data type with minimal loss of information. The net rate of process- 
ing scalar measurements is then Mf, where M is the number of data types and f is 
the filter update frequency. The possibilities for minimizing M are limited, and 
a value of 1 Hz was found satisfactory for f in the simulation tests described 
next. 


The influence of filter update interval on accuracy was of considerable interest 
in the prospective flight application. Between updates, the filter is an unaided 
inertial navigator with error divergence behavior which reflects its acceleration 
estimation accuracy. This accuracy is sufficiently low in the present system during 
maneuvering that performance losses become significant for update intervals in the 
range of 1 to 10 sec. The nature of these effects is illustrated by the comparison 
of sample case histories in figure 3.2 for intervals of 0.5 and 10 sec. Position 
and velocity error histories are sawtooth-like functions which drift during the update 
interval relative to the error histories corresponding to continuous or high-rate 
measurement processing; they drift according to the equations 

a dt 2 


t < t < t + A 
o o 

'“O 

and are approximately returned to the continuous processing error history discretely 
at a filter update. The amount of drift depends on the update interval, A, the 
velocity error just after the previous update, V D , and the acceleration error, a. 

In figure 3.2(a) the lateral axis results show significant position and velocity error 
drifts for the 10-sec interval; drifts for the normal axis are much smaller (after the 
initial transient) as a result of much better acceleration accuracy. In addition, the 
lag and loss of information at larger intervals results in poorer transient response, 
as is evident in the degraded initial transients for all states and the significant 
loss of lateral acceleration accuracy during the turn (legs 4 and 5), where maneuver- 
induced gyroscopic error transients occur. The performance effects of practical 
interest are those on trajectory tracking errors and control activity. A generic, 
automatic, trajectory tracking system was included in the simulation, and its response 
to the sample case estimation errors is shown in figure 3.2(b). Several effects are 
visible, particularly for the lateral axis. First, estimation error jumps at each 
update result in corresponding control command histories (calculated as a corrective 
acceleration command) that are increasingly characterized by rate and authority limit 
saturation as the interval increases, and, second, tracking errors show increased 
excursion extremes and lags. Thus, both tracking accuracy and control activity 
degrade significantly with increasing interval over the range of intervals tested 
here. 


t Q +A 


R(t + A) = R + V A + 
o o o 


// 


t 0 +A 

V(t + A) = V + I a dt 
o o J 


Finally, Monte Carlo simulation results for some effects of filter update inter- 
val are shown in figure 3.3. These data are taken from a turning segment of the STOL 
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A = 0.5 sec 



(a) ESTIMATION ERRORS (b) TRAJECTORY TRACKING ERRORS AND 

CONTROL ACTIVITY 

Figure 3.2.- Effects of update interval on estimation and control performance: 

sample case, STOL approach trajectory. 

approach where acceleration estimation accuracy is poorest. Root-mean-square jumps 
in the estimate at the update increase nearly linearly with update interval, and the 
corresponding effects on tracking accuracy and control activity (not shown) are also 
proportional to the interval. Root-mean-square estimation errors at the end of the 
interval indicate the poorest accuracy that occurs during a filter cycle; these errors 
are seen to be insensitive to the interval size below intervals of 2 to 3 sec. This 
occurs because some error drifts can be temporarily in the direction of smaller mag- 
nitudes, but this effect is lost as interval size increases. Thus, these data indi- 
cate that intervals of 1 sec can be selected with no significant increase in the 
estimate jumps or loss of estimation accuracy compared with measurement processing at 
higher rates, and this choice has been found satisfactory generally. Larger intervals 
may also be satisfactory, but further increases in interval size are decreasingly 
effective in reducing the required computation time. 

The algorithm outlined above was simulated on a CDC 7600 digital computer, and 
an equation summary for this simulation is provided in table 3.7. These equations 
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UPDATE INTERVAL, sec UPDATE INTERVAL, sec 

(a) ESTIMATE JUMP SIZE AT UPDATE (b) ACCURACY AT END OF UPDATE INTERVAL 

Figure 3.3.- Effect of filter update interval on estimation accuracy: 

STOL approach, turn segment. 


repeat those previously given in this section and require little additional comment. 
We note that <J>, {Vf^}, Q are sparse matrices and that only their nonnegligible 
elements are stored in the implemented algorithm along with arrays of indices which 
permit the elimination of trivial multiples. The simulation includes the computation 
lags associated with real time operation (0.4 sec for the filter equations). Further 
details on these lags and on coordination of the multirate computations in real time 
are given in reference 26. We also note that other square-root covariance algorithms 
have been proposed in the literature since Potter’s algorithm was selected for this 
work and routine developed. Some of these offer potential reductions in computation 
time and their evaluation for real-time, fixed or floating computations would be of 
interest. However, navigation accuracy would not be affected since it is insensitive 
to measurement reception and processing rates at the rates used in the present 
algorithm. 
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Summary 


The existing theory and applications devices have been applied to obtain a 
Kalman filter algorithm appropriate to the context of this study. The principal 
design issue is to realize the maximum available estimation accuracy from the set of 
sensors that is used while minimizing the required computation time in a real time 
flight computer implementation. The required computation time depends principally on 
the rate of processing scalar measurements with the Kalman filter; this rate is mini- 
mized with negligible loss of information by accumulating measurements over an inter- 
val and compressing them to a single, nearly equivalent, scalar measurement for each 
data type which can then be processed by the filter once per interval. For this 
purpose it was found empirically that measurements could be sampled and accumulated 
at 10 Hz with negligible loss of information and processed at 1-sec intervals without 
significant loss of estimation accuracy. Measurement sample error correlation times 
and the IMU accuracy are important factors in establishing these values. This device 
provides an order of magnitude reduction in the rate of processing scalar measure- 
ments from that which would be required without measurement compression. 

Additional design factors affecting computation time and estimation accuracy 
remain to be studied in simulation tests described in the following sections. These 
include (1) evaluation of several candidate separations of the system and filter into 
independent lower-order systems and filters to obtain significant reductions in the 
time required to process a scalar measurement, and (2) appropriate choice of acceler- 
ation measurement error states from the three candidate formulations noted in this 
section in order to improve accuracy. The selection of state variables is otherwise 
largely dictated by control requirements and by a brief review of the parameters of 
the measurement functions, their a priori accuracy, and the availability of signifi- 
cant information from the sensors compared with the a priori information. 


4. ATTITUDE AND ACCELERATION ESTIMATION ACCURACY 


The object of this section is to determine the best of several candidate formu- 
lations of the acceleration-measurement error states, the acceleration and attitude 
accuracies achieved by the filter, and the sensitivity of the acceleration and atti- 
tude accuracies to sensor accuracy. 

Three formulations of the acceleration measurement errors will be studied; these 
are given in terms of (I) runway axes measurement biases (14-state filter), (2) atti- 
tude measurement biases and a vertical acceleration measurement bias (15-state fil- 
ter), and (3) level heading axes measurement biases (modified 14-state filter). 

These formulations are approximate models of the low frequency error processes 
that result in acceleration measurement errors; the principal error sources are the 
fixed accelerometer errors (biases, misalignments, and scale factor errors) and 
maneuver-induced error transients in the attitude gyroscopes. The 15-state filter 
models the gyroscope error sources directly and is of interest for possible accuracy 
increases in attitude, as well as acceleration estimation. The modified 14-state 
filter uses an axis frame in which the acceleration measurement errors due to the 
fixed accelerometer errors are nearly constant, thus removing one source of the time 
variations found in the runway axis error states. 
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Reference Trajectory and Ensemble of Approaches 

Accuracy results are given as rms estimation error histories for a 10-sample 
ensemble of approaches along a reference trajectory. The reference trajectory is 
defined in figure 2.1 and table 2.1 from terminal area entry to landing. It contains 
the range of maneuvers, accelerations, speeds, and flightpath angles expected in STOL 
approach operations, including a speed range from 140 to 65 knots at landing, and a 
decelerating, helical descent and glide slope at -7.5°. The path is coarsely defined 
as a sequence of straight-line and circular-arc legs for which kinematic data are 
listed in table 2.1. The coarsely defined trajectory is discontinuous in velocity 
and acceleration at the leg junctions, but the trajectory generating algorithm (given 
in ref. 27) smooths this to a continuous trajectory which satisfies appropriate oper- 
ational and aircraft limits on velocity excursions, acceleration, and jerk. 

Measurement biases, noise, and other error types for the ensemble of 10 approaches 
are sampled randomly or generated in accordance with the statistical and deterministic 
models previously defined. The wind environment is constructed by superposition of 
mean and turbulent wind components generated by models commonly used in aircraft simu- 
lations for the low-altitude winds. The direction and magnitude of the mean wind are 
selected randomly for each approach such that its 3 a values match those recommended 
by the FAA for the landing zone as a function of direction (25-knot headwind, 15-knot 
crosswind, 10-knot tailwind). The mean wind is assumed constant for periods much 
longer than the approach flight duration and increases with altitude by a factor of 2 
in the atmospheric boundary layer (up to 300 m) and is constant above that. The wind 
turbulence history is that encountered by an aircraft flying through a space-frozen 
turbulence field and is generated as a three-dimensional Gaussian random process 
(Dryden model) with zero means and variances that depend on aircraft axis, attitude, 
and airspeed, and that are proportional to the field intensity (field intensity is 
measured as the landing zone horizontal gust variance and is a random variable with a 
mean value of 0.7 m/sec). Documentation of this wind model is omitted here, but an 
extensive literature on this subject exists, including theory, observations, and 
recently recommended models for terminal area aircraft simulations (e.g., refs. 28 
and 29). The aircraft motion simulation omits aircraft attitude dynamics and turbu- 
lence response so that their effects on actual aircraft accelerations and some higher 
frequency sideslip and attitude angle variations are absent; their effects on navi- 
gational accuracy are noted where significant. 

The randomly selected constants for each of the 10 sample approaches are listed 
in table 4.1; these sample values are used in all ensemble results of the report 
unless otherwise stated. In this section, interest focuses on performance during the 
maneuvering segments of the trajectory. Empirical results are, therefore, given 
beginning near waypoint 2 (t = 647) where maneuvering starts. 


Attitude Estimation 

Attitude estimation accuracy for the 15-state filter is compared with the accu- 
racy of the gyroscope, $g, 0g, $ g, in figure 4.1. A legend indicates the trajectory 
leg number and navaids in use. These results show no success in improving attitude 
accuracy. Roll and pitch estimation accuracies are considerably poorer than the 
errors in the vertical gyroscope, and the heading error is nearly identical to the 
error in the directional gyroscope; that is generally. 
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rms(<j>) > rms($ ) 
6 

rms(0) > rms(0 ) 
S 

rms(ij)) = rms(i(j ) 
§ 


This occurs in part because accelerometer errors are significant and must be inter- 
preted by the 15-state filter as equivalent attitude measurement biases. Both equa- 
tion (2.14) and the empirical results for sample cases show that longitudinal and 
lateral accelerometer biases are interpreted, respectively, as pitch and roll measure- 
ment biases. In addition, equation (2.14) shows that the principal deterministic 
effect of \pg on am is longitudinal, nonzero only in turns, and indistinguishable 
from the effects of 0g, which have a much larger gradient throughout the present 
test trajectory. The result is that the filter interprets the effect of ij/g almost 
entirely as an equivalent pitch bias; this is evidenced in figure 4.1, both by the 


TABLE 4.1- RANDOMLY SELECTED CONSTANTS FOR ENSEMBLE OF APPROACH FLIGHTS 


Sample 

Wind parameters 

Measurement biases j 

Mean wind , 
m/sec 

Turbulence 

intensity, 

m/sec 

V0RTAC 

Baro- 

M0DILS 

Body-mounted 
accelerometers, g 

Range, 

m 

Bearing, 

deg 

altimeter, 

m 

Range, 

m 

Azimuth, 

deg 

Elevation, 

deg 

Longi- 

tudinal 

Lateral 

Normal 

w x (0) 

y°> 

1 

-1.2 

-7.6 

mem 

252 

0 

-11 

-3.7 

0.02 

-0.046 

0.012 

-0.004 

-0.042 

2 

4.0 

6.6 


-121 

-1.0 

17 

10.5 

.15 

.016 

.017 

.010 


3 

.9 

-7.1 


282 

-2.5 

41 

-9.9 

-. 15 

.015 

.014 

-.031 

-.018 

4 

9.2 

7.3 

1.0 

-134 

-2.9 

-15 

4.4 

-.11 

-.027 

.018 


.007 

5 

9.8 

1.8 

.7 

61 

-1.3 

4 

1.5 

.11 

-.008 

.011 

.032 


6 

13.0 

-2.1 

.8 

-93 

-.4 

24 


.17 

-.031 

-.003 



7 

2.7 

5.9 

.5 

-348 

-.6 

67 

-12.8 

0 

.016 

-.016 


.004 

8 

.5 

3.6 

1.3 

-386 

-2.4 

-14 

-2.0 

.14 

.041 


-.025 

-.004 

9 

5.5 

.4 

1.0 

2 

-.3 

-27 

4.1 

.12 

.055 



.004 

10 

14.6 

-3.7 

1.4 

233 

2.6 

6 

1.9 

.23 

.051 

mSm 

iSS 

-.010 

Mean 

5.9 

0.5 

1.0 

-76 

-0.9 

9 

-0.9 

0.07 

0.008 

0.003 

-0.005 

-0.006 

RMS 

7.9 

5.2 

1.1 

226 

1.7 

32 

6.6 

0.15 

0.035 

0.0125 

0.018 

0.016 
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Figure 4.1.- Attitude estimation accuracy of 15-state filter (rms errors). 

behavior of ifi and the (delayed) correspondence of some variations in 0 with 
variations in $ g . 

Additional results with errorless accelerometers were obtained to determine if 
improvements in this sensor would permit the filter to estimate gyroscopic errors to 
a useful degree. The results (fig. 4.1) show no success; roll angle is unambiguously 
detectable in this case but is improved only to about the accuracy of the gyroscopic 
roll angle, and the pitch error excursions, which are well above 0g, owing to misin- 
terpretation of the effects of $g, remain present, and is virtually unchanged. 

The filter’s information on attitude errors is derived from position measure- 
ments; figure 4.1 shows similar performance using either VORTAC or MODILS and suggests 
that performance is insensitive to position measurement accuracy in the accuracy range 
of these navaids. This is confirmed by results using errorless position measurements 
(not shown); in this case, roll estimation accuracy is significantly improved, with 
peak excursions to 0.7°, but pitch and heading remain subject to the error effects 
noted above and are unimproved. 

These results show that the present sensors, in combination with an appropriate 
formulation of the filter states, cannot improve the accuracy of the measurements of 
the attitude gyroscope. 


44 




Acceleration Estimation 


The estimation accuracy for the path axes coordinates of acceleration is shown 
in figure 4.2, where it is compared with the accuracy of the measured acceleration 
derived from instantaneous sensor outputs. As seen, the two filters perform simi- 
larly. Both succeed in reducing measurement errors significantly as a result of 
their ability to estimate and remove some effects of the fixed and low frequency 
accelerometer and gyroscopic errors; both give good, nearly invariant normal axis 
accuracy in excess of 0.01 g after the initial transient; and both give similar 
lateral axis accuracy. For the longitudinal axis, the error reduction is modest and 
the performance of the 15-state filter is poorer, with several excursions in excess 
of the measurement error that can be traced to excursions in 0g,$g. 

Accuracy for the vertical axis and horizontal plane components of acceleration 
differ significantly because (1) the first-order effects of gyroscopic error tran- 
sients are confined to the horizontal plane (eq. (2.14)) and (2) the fixed acceler- 
ometer errors are nearly constant along the vertical axis because of the low roll and 
pitch angles used in passenger operations; however, the latter are time-varying rela- 
tive to the horizontal plane runway axes during turns. Thus, the vertical accelera- 
tion measurement error is nearly constant and the filter succeeds quite well in esti- 
mating this from position measurements. In contrast, the horizontal plane components 
are both larger and time varying during maneuvering flight and this continually 
excites the filter's transient response. It is the variation rate of the measurement 
errors rather than their magnitude that limits the accuracy obtained here. 

In view of the results seen in figures 4.1 and 4.2, the 15-state filter can be 
dropped from the discussion; the increased number of states adds a poorly observable 
state, ip , to the computations, without increasing the accuracy of the estimated 
acceleration . 
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Figure 4.2.- Acceleration estimation accuracy: 14- and 15-state filters. 
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The separate influence of the attitude gyroscope and accelerometer errors is 
examined in figure 4.3, in which the existing acceleration estimation accuracy is 
compared with that achieved with (1) errorless gyroscopes and (2) errorless acceler- 
ometers. In each case, the magnitude of the estimation error indicates the contri- 
bution of the normal sensor type to the total error, as well as the maximum perfor- 
mance available from sensor improvements. For the longitudinal and lateral 
components, estimation errors owing to accelerometer errors are close to the total 
error over most of the trajectory, and those owing to gyroscopic errors are usually 
somewhat smaller; performance for both components is poorer than for the normal 
axis and shows excursions above 0.01 g over some or most of the trajectory. Thus, 
accelerometer errors dominate the estimation error for this filter. 

The same comparison is made in figure 4.3(b) for the modified 14-state filter 
which uses level heading axes acceleration measurement bias states. As seen, acceler- 
ation estimation errors for the longitudinal and lateral axes appear to be almost 
entirely the result of gyroscopic errors (except initially); the variations in estima- 
tion errors follow the gyroscopic error histories of figure 4.1 in accordance with 
equation (2.14). This result is quite different from that obtained in figure 4.3(a) 
and indicates much poorer performance for the modified 14-state filter in detecting 
the effects of gyroscopic errors. The cause of this difference in performance is not 
clear. On the other hand, the modified 14-state filter shows excellent performance 
in estimating accelerometer biases when the gyroscopes are errorless (accuracy is 
invariant and exceeds 0.01 g on all axes), as was anticipated. 

Comparing the two 14-state filters, the net result is that each filter performs 
better with one type of sensor error but does poorly with the other. The accuracy 
achieved by the modified 14-state filter is nearly identical to that of the 15-state 
filter, which also has the potentially favorable property that the^fixed accelerometer 
errors map into nearly constant state errors (in this case, <j>, 0, 2). Thus, a com- 
parison of the two filters would give the same result as that shown in figure 4.2, 
and we select the runway-axis formulation of the acceleration bias states for the 
filter as a result of its somewhat better accuracy, which is due to its greater (but 
limited) ability to detect the error effects of the attitude gyroscopes. 

The effect of position measurement accuracy on acceleration estimation accuracy 
was also examined. In figure 4.2 both VORTAC and MODILS data result in the same range 
of minimum to maximum errors during maneuvering So that performance appears insensi- 
tive to position measurement accuracy over the range obtained from these navaids. 
Further tests were made comparing performance for errorless and normal position mea- 
surements. First, results were obtained for a straight line static equilibrium path, 
using VORTAC data and large accelerometer biases (0.05 g) ; the ensemble averages 
showed that VORTAC gave good steady state, acceleration estimation performance (below 
0.01 g) — equivalent to that for the normal axis in figure 4.2. Thus, under the 
restricted circumstances of inertially constant acceleration measurement biases, even 
fairly poor position measurements suffice to calibrate the present IMU to better than 
0.01 g. The errorless position data, therefore, yielded little improvement in steady 
state accuracy, but significantly reduced the filter’s settling time in estimating a r . 
Secondly, the same comparison was made for. the STOL approach path to determine if the 
faster response time would maintain the acceleration estimation accuracy closer to 
the steady state performance during maneuvering. The results (fig. 4.4) show moder- 
ate improvements throughout the approach in longitudinal and lateral acceleration, 
with rms errors of 0.02 g or better, but accuracy remains time-varying in turns. 
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Figure 4.4- Acceleration accuracy: effect of position-sensor accuracy. 

Other simulation tests (not shown) indicated that acceleration estimation per- 
formance of the three filters was insensitive to the airspeed sensor accuracy and to 
the values of the correlation time-constants used in the filter's state equations for 
acceleration-measurement errors (table 3.2). These models are tractable but fic- 
tional, and their parameters are "tuned" to achieve the best performance with the 
simulated sensor error models. Although it is best to tune a flight system with 
recorded flight data or with well-verified simulation models, the observed insensi- 
tivity to these time constants implies a related insensitivity of the in-flight 
accuracy to simulation model inaccuracies. 


Summary 

The simulation tests discussed in this section lead to six major conclusions 
regarding attitude and acceleration estimation. 

1. Three formulations of the acceleration measurement error states were tested 
and the formulation as runway axis bias errors was selected, based on its somewhat 
better acceleration estimation accuracy on a maneuvering STOL approach. 

2. Attitude estimates can be given only to the accuracy of the attitude gyro- 
scopes. Efforts to improve on this by formulating the filter states as gyroscope 
errors resulted in even poorer accuracy because (1) the fixed accelerometer errors 
are interpreted as equivalent gyroscope errors, and (2) 0g and $g are indistinguish- 
able as causes of position residuals. Further, tests with errorless accelerometer 
and position data showed that any success deriving from these sensor improvements 
was limited to roll angle. 

3. The dynamics of the attitude gyroscope were continually excited by maneuver- 
ing on the final part of the STOL approach path until the glide slope segment; this 
resulted in rms pitch and roll errors between 0.5° and 1.5° and rms heading errors 
between 1.5° and 5° during most of the time. 
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4. The filter succeeds in estimating acceleration more accurately than acceler- 
ation can be measured by using position data to detect the fixed and low frequency 
acceleration measurement errors. Accuracy for the vertical axis is good, with nearly 
stationary estimation errors below 0.01 g rms throughout the approach; this is a 
result of the easily detected, near-constant acceleration measurement bias for this 
axis. Accuracy for the horizontal plane components is significantly poorer, with 
time-varying estimation errors resulting from the dynamics of the attitude gyroscopes 
and the time-varying runway axis components of the fixed accelerometer package errors 
during maneuvering. The accuracy range on the test approach path was as follows: 


Component 


Estimation error, rms 


Longitudinal 

Lateral 

Normal 


0.01 to 0.03 g 
0.005 to 0.025 g 

<0.01 g 


5. Tests showed that the fixed accelerometer package errors were the principal 
source of the horizontal plane acceleration estimation errors. An alternative formu- 
lation of the acceleration bias error states in level heading axes succeeded in 
detecting the fixed accelerometer errors to better than 0.01 g but was unsuccessful 
in improving accuracy because of reduced performance in detecting gyroscope errors. 


6. Significant improvements over the acceleration and attitude accuracies 
achieved here do not appear possible by further development of the filter formulation 
nor from improved accelerometer or position sensors; such improvements can be 
obtained, however, from a more accurate IMU, such as an inertial platform or Schuler- 
tuned strapdown system (ref. 30) or, perhaps, by adding rate gyroscopes to the present 
set of sensors. 


5. POSITION, VELOCITY, AND WIND ESTIMATION ON A STRAIGHT LINE PATH 


Performance details are studied next, 
using results from a level, straight 
flightpath near the runway. In these 
tests, maneuver-induced acceleration mea- 
surement error transients are absent and 
acceleration accuracy exceeds 0.01 g on all 
axes. This simplifies somewhat the analy- 
sis of the dependence of estimation accu- 
racy on sensor accuracy and flightpath. 

The flightpath (see sketch E) is parallel 
to the runway at 2,200 m and is flown at 
110 knots; the filter is initialized out- 
side MODILS coverage and acquires MODILS 
after 200 sec. 

The objectives of this section are to 
determine estimation accuracy, to determine 
its principal transients on unaccelerated 
paths, to determine its sensitivity to 
sensor accuracy, and to settle empirically 
several filter design issues. Performance 


x ~ km 
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details are reviewed separately for each type of state and for states associated with 
the horizontal plane and vertical axis motions. Ensemble rms errors for all states 
of interest here are given in figure 5.1. Note that scales are changed at the switch 
to MODILS data for some states in order to improve resolution. 



0 100 200 300 

TIME, sec 


VORTAC 

MODILS 

BARO 




STANDARD CASE 

HIGH ACCURACY IMU (a a = 0.001 g) 

WIND STATES DELETED 


Figure 5.1.- Performance comparisons for (1) standard case, (2) high-accuracy IMU, 

(3) filter without wind states. 


V0RTAC Biases 

The V0RTAC bias states are estimated as zero at filter initialization so that 
their initial sample and ensemble ms errors equal the sample case and average biases, 
respectively. Subsequently, information on these biases is negligible on the straight 
line path compared with the a priori information, until (1) passage close to the VORTAC 
site (during 100 sec < t < 150 sec in fig. 5.1) permits calibration of the range bias 
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to an accuracy of 60 m from acceleration measurements; and (2) entry to MODILS cover- 
age (at t = 200 sec) permits calibration of both VORTAC biases to the accuracy of the 
MODILS measurements (15 m, 0.25°). After these events, the VORTAC remains calibrated 
to the accuracy of MODILS regardless of later switches between MODILS and VORTAC so 
that no significant maneuvering to an offset trajectory is excited by the estimator 
at such switches. 

Before receiving MODILS data, information on the bias states can be obtained 
theoretically from (1) acceleration measurements, and (2) mutual calibration of the 
range and bearing to, at best, the accuracy of the better distance measurement. These 
possibilities are explained next, but the empirical results show only limited success 
from these sources. 

First, suppose that position is calculated from biased VORTAC measurements. The 
acceleration required to follow this calculated position differs from the accelera- 
tion that is measured along the actual path and the difference provides nonnegligible 
information on the VORTAC biases, provided the acceleration measurements are suffi- 
ciently accurate. More explicitly, the error in computing position from noiseless 
VORTAC measurements with unknown biases, b tr ,b tb is 

R (b) = b tr u + r t b tb £ (5.1) 

where (jJjj>) are the unit radial vectors 
from the VORTAC and its perpendicular (see 
sketch F). If, now, velocity and acceler- 
ation are calculated as derivatives of the 
computed position, the errors that are due 
to the unknown bias are given by deriva- 
tives of equation (5.1): 

-(b) - 

- = b tr\£ + b tb^r ® 1 (5.2) 

' = b tr ($ t £. ~ 'f'^u) + b tb k r ® a (5.3) 

where is the bearing angle. Bearing 

bias error results in (1) an apparent 
velocity that is lateral to the actual 
aircraft velocity and (2) an apparent 
acceleration error that is lateral to the 
actual aircraft acceleration. The apparent 
acceleration b^a is zero for the present 

unaccelerated straight flightpath so that bearing bias error is unobservable to the 
acceleration measurements here. The estimator, however, will settle on parameter 
values such that the VORTAC measurement residuals have zero mean; this can be any 
combination of errors such that 



Sketch F 


E[V+V (b) ]=0 (5.4) 

Using equation (5.2), the result for the lateral axis is 

E [V * Ip + b tb V + b tr ^ t £ • i p ] = 0 (5.5) 
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In most situations, the term in b tr is negligible, so that the bearing-bias error 
results in a corresponding lateral velocity estimation error bias that satisfies 

E[V * Ip] = “V E[b tb ] (5.6) 

The initial errors (V(0), btb(O)) are independent and do not satisfy equation (5.6) 
in general, but their steady state mean values are driven to this line. This can 
result in noticeable sample case transient changes in the bearing bias estimate, even 
though the ensemble average for b t b changes very little in figure 5.1.^ This is 
demonstrated in figure 5.2(a), which shows the trajectories of { ( V • ^pjbtb)^ ^ 0} 
from selected initial errors. In each case, the figure gives the mean trajectory for 
an ensemble of samples having the restricted initial errors. As seen, the combined 
bearing bias and lateral velocity estimation errors are driven to the line defined by 
equation (5.6). 

On the other hand, if we consider accelerating flightpaths, then the apparent 
acceleration |b t ^a| is largest during turns when it is oriented along the level 
longitudinal axis (k r ® a. = g tan (J> iQ , but is, at most, only 0.006 g per degree of 
bearing bias in passenger operations (<f) < 20°). This provides only crude bearing cali- 
bration accuracy of the same size as its a priori accuracy (2°) from an IMU 
with an accuracy of 0.013 g, but a substantially better calibration can be achieved 
by an inertial grade IMU. Figure 5.2(b) shows results for both the standard and high 
accuracy IMU during a turn at a distance of about 23 n. mi. from the VORTAC station. 

As seen, the standard IMU achieves some improvement in spite of poor acceleration 


E[V • jp] m/sec 



RESIDUALS: 


E[b tb l = 



(a) Behavior of bearing bias and lateral velocity estimation errors. 
Figure 5.2.- VORTAC bearing-bias estimation outside MODILS coverage. 
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STRAIGHT | TURN STRAIGHT 


SEGMENT 

(b) Bearing calibration by IMU during a turn (V = 150 knots, R c = 1675 in), 


Figure 5.2.- Concluded. 


estimation accuracy, and the inertial grade IMU detects the bearing bias to an accu- 
racy better than 0.5 ; this results in useful reductions of position errors and of 
the lateral velocity bias error due to bearing bias error. Results with the inertial 
grade IMU on a straight line path at the same distance are included in the figure and 
demonstrate that the turn is required, as well as the inertial grade IMU, to obtain 
these performance benefits. Note that the required turn can be executed early in the 
approach, since the available calibration accuracy is independent of distance to the 
VORTAC. However, the capacity to calibrate angle measurements with this IMU is 
limited and does not suffice to similarly improve the a priori accuracy of the MODILS 
azimuth calibration. 


The apparent acceleration due to range bias error is negligible in equa- 
tion (5.3), except during periods of nonzero bearing rate, ip t ; for rates of 1.75° /sec, 
which occur when passing close to the VORTAC station in the present example, the rms 
acceleration error resulting from range bias is 0.03 g along the VORTAC radial. This 
permits calibration to 100 m, with 0.01-g acceleration estimation accuracy and 
accounts for the improved range calibration accuracy seen in figure 5.1 before enter- 
ing MODILS coverage. More generally, higher bearing rates result in greater calibra- 
tion accuracy, but this effective opportunity to calibrate the range measurement 
occurs only when passing close by the VORTAC station. 

The range and bearing bias errors can be functionally related by taking VORTAC 
measurements at different times along any path that is not a VORTAC radial. If one of 
the measurement types provides a sufficiently more accurate distance measurement than 
the other, the estimated bias of the poorer measurement can be improved. For example, 
near the VORTAC station, bearing provides a distance measurement 3 or 4 times more 
accurate than the a priori range bias estimate, but at a distance of 20 n. mi. the 
converse is true. However, this source gave only negligible information on the bias 
errors compared with the a priori or IMU information and we omit an analysis. 


53 




Horizontal Plane Position Coordinates 


o 


The initial estimate is computed from VORTAC data so that initial errors are dic- 
tated by VORTAC measurement errors and depend on the location where filter initializa- 
tion occurs. In the present test, initialization occurs close to the VORTAC site and 
rms errors are of the order of 250 m in both directions. At other locations, initial 
errors would be unchanged in the direction of the VORTAC radial, but would increase 
with range for the direction lateral to this and are much larger (0.75 n. mi.) at 
entry to the terminal area. After initialization (fig. 5.1), accuracy changes very 
little until passage close to the VORTAC permits calibration of the range bias to 60 m 

and a corresponding reduction in rms 
position error below 100 m in both 
coordinates; this is followed by a 
rapid reduction by an order of magni- 
tude at entry to M0DILS coverage. 

Before these events, the position esti- 
mation error is biased by the VORTAC 
calibration error vector (eq. (5.1)) 
and sample case position errors are 
larger or smaller, depending on the 
sample case VORTAC biases. 

Sensitivity to acceleration mea- 
surement accuracy is shown in fig- 
ure 5.1 by the comparison with results 
for a high accuracy IMU (c a = 0.001 g) . 
Improvements for the position states 
are limited to an improved capacity to 
detect range bias when passing close 
to the VORTAC and related local 
improvements in position accuracy. 



iRm I 


m 


(a) POSITION ESTIMATION ACCURACY 



POSITION MEASUREMENT ACCURACY - m 


The dependence of position estima- 
tion accuracy on the combined accuracy 
of both position and acceleration mea- 
surements is examined further in fig- 
ure 5.3(a). These results were 
obtained in a series of simulation 
tests; the estimation accuracy was 
obtained from the covariance matrix as 

cj R | = E[(x 2 + y 2 )] = a 2 + a 2 

and measurement accuracy was varied 
over the range of interest, using 
unbiased MODILS measurements with error 
and flightpath parameter values 
selected to give equal accuracy in all 
horizontal directions in the expression 


(b) VELOCITY ESTIMATION ACCURACY 



= E[(Y u + 
L v mr— 


¥ p) 2 ] 

xym ma 1 - 7 


Figure 5.3.- Effect of position and acceler- = a 2 + d 2 a 2 

ation measurement accuracy on performance. Tnr x Y m ma 
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The results indicate several conclusions. First, estimation accuracy is substantially 
better than the position measurement accuracy in all conditions — a result of the 
filter’s smoothing of measurement noise. Second, improvements in either acceleration 
or position sensing accuracy yield some reduction in estimation errors, but greater 
reduction is achieved by reducing position sensing errors than by reducing accelera- 
tion sensing errors by the same factor. Two points, labeled "standard" case and 
"high-accuracy IMU," indicate that only a moderate improvement in position estimation 
accuracy would be realized by replacing the present accelerometer attitude-gyroscope 
system with an inertial grade IMU. Thus, figures 5.1 and 5.3(a) both indicate that 
position estimation accuracy depends principally on the position measurement accuracy, 
and significant improvements in position accuracy over the present performance 
requires improved position sensors. 


Horizontal Plane Velocity Coordinates 

The velocity estimate is initialized at the measured air velocity vector, (Y^ a , 
Yya» 0) , and its initial error is dominated by the actual wind with a theoretical 
standard deviation of 6.1 m/sec in each coordinate. Subsequently (fig. 5.1), the 
velocity accuracy settles in about 25 sec to about 2 m/sec, shows a modest transient 
when passing close to the VORTAC, and settles to 1 to 1.5 m/sec after acquiring 
MODILS data. For comparison, the airspeed sensor accuracy is 0.6 m/sec, and this 
sensor has been used in restricted control modes (e.g., automatic airspeed capture 
and hold modes and airspeed-referenced trajectory tracking). Accuracy for the hori- 
zontal plane inertial velocity coordinates with the present sensors and estimator is 
significantly poorer than this, except under conditions of maximum position and 
acceleration measurement accuracy. 

More generally, the transient and steady state behaviors of the velocity accuracy 
show several features. First, the initial settling time depends somewhat on position 
accuracy and direction; it is faster using MODILS data, slower at large distances for 
the direction lateral to the VORTAC radial, and is faster for the longitudinal axis 
than for the lateral axis because of aiding by airspeed measurements. 

Second, the steady state lateral velocity error, using VORTAC data, is moderately 
biased by the bearing bias error. Bearing bias is very poorly observable before 
entering^MODILS coverage, but its presence is interpreted as a lateral velocity in the 
amount b^V in accordance with equation (5.2). This effect is illustrated in 
figure 5.4(a), which shows the ensemble means for velocity in the restricted case in 
which the bearing is 4° in all samples; as seen, the lateral velocity error, y, has a 
steady state bias (4 m/sec) while the longitudinal velocity error, x, is unbiased. 

The steady state velocity accuracy is otherwise independent of the remaining, ran- 
domly sampled constants of each flight (mean winds and measurement biases). 

Third, significant error transients when passing close to the VORTAC site can 
occur in some samples. Equation (5.2) indicates that estimation errors in VORTAC 
range bias, velocity lateral to the VORTAC radial, and acceleration along the radial 
are all possible causes of range measurement residuals. If one of these errors is 
significantly larger than its standard deviation in the filter’s covariance matrix, 
then the filter temporarily misinterprets the residuals when near the VORTAC and pro- 
duces error excursions in all of these states. This is illustrated in figure 5.4(b) 
by the ensemble means in the case in which range bias is restricted to its 2a value 
(610 m) in all samples; as seen, velocity accuracy shows a significant transient. 

The remaining source of significant velocity accuracy transients is the IMU errors, 
and they are discussed in a later section. 
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(3) (t» tr , b tb> = <°. 4 °) (b) (b tr ,b tb ) = (610m,0) 


Figure 5.4.- Velocity and wind estimation: effects of large VORTAC calibration 

errors. 

The dependence of velocity accuracy on position and acceleration measurement 
accuracies is shown in figure 5.3(b). Velocity estimation accuracy is obtained from 
the filter’s covariance, using 


o 


2 



2 , 2 
a . + a . 
x y 


and the results correspond to white Gaussian sensor errors. These results show sig- 
nificant dependence on the accuracy of both types of measurements. Two points, 
labeled "standard" and "high accuracy IMU" are noted in the figure. They indicate 
that a substantial improvement, to 0.6 m/sec, would be achieved with MODILS data by 
using an inertial grade IMU. More generally, to achieve a velocity estimation accu- 
racy of 1 m/sec throughout the terminal area requires improvements in both the 
position and acceleration measurement accuracy of the present system. Improved posi- 
tion sensing is needed to remove the transient and steady state effects of the VORTAC 
biases discussed above, and a high accuracy IMU is needed to eliminate the effects of 
acceleration measurement error excursions excited during maneuvering segments of the 
approach. A high accuracy IMU by itself cannot achieve this performance for lack of 
the ability to estimate the VORTAC biases accurately (as shown in fig. 5.1), and an 
improved position sensor by itself lacks the ability to estimate acceleration mea- 
surement errors on maneuvering paths (as shown in the previous section) . 


Wind States, Airspeed Measurements, and Sideslip Effects 

Because the wind is initially estimated as null, initial errors equal the actual 
winds and have a standard deviation of 6.1 m/sec. Subsequently, accuracy (fig. 5.1) 
is nearly identical to that of the velocity coordinates. This gross behavior can be 
understood from the approximate wind estimate, W' , obtained from the current velocity 
estimate and derived air velocity measurement. 
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with estimation errors that have zero mean and variances 



a? 

x 


+ 




a? + 

y 



(5.7) 


In most of the present test, these errors are dominated by the inertial velocity 
terms, which range from 6.1 m/sec initially to 1.5 m/sec in steady state on MODILS 
data; this compares with the theoretical air velocity measurement accuracy of 
1.8 m/sec. Thus, wind estimation accuracy and its transient and steady state response 
to initial errors and VORTAC biases, and its sensitivity to position and acceleration 
measurement accuracy are largely the same as for the velocity coordinates previously 
discussed. In cases in which the velocity estimate is much more accurate than the 
air velocity measurement, as in the high accuracy IMU case with MODILS data shown in 
figure 5.1, the air velocity measurement accuracy limits the steady state wind esti- 
mation accuracy. 


The treatment of wind estimation was an issue in the filter design. The inclu- 
sion of wind states in the filter permits processing air velocity measurements, which 
improves output accuracy both normally and during dead reckoning in the event of 
position navaid dropout. Alternatively, a substantial saving in computation time can 
be obtained by removing the wind states and air velocity measurements to a separate 
filter, as shown in sketch G. This can bi 
done without affecting the observability 
of the remaining states. The effect on 
accuracy is examined in figure 5.1, which 
includes results for both the normal and 
separated filters; their comparison shows 
significant accuracy losses in all states 
associated with the horizontal plane 
motion, principally during VORTAC use. 

Thus, the separated filter discards sub- 
stantial information on these states; 
because it does, we retain the wind state; 
in an integrated filter. Sketch G 

Crosswind is indistinguishable from sideslip angle and heading error for the 
present sensors. The air velocity measurement is the only measurement with functional 
dependence on the wind. Recalling equations (3.13) and (3.16) (which are repeated 
here) , the air velocity measurement is derived from airspeed and heading measurements 
as 




Va 


cos ^ g \ 

sin v 


and represented in the filter as 






(3.13) 
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with Gaussian white errors. However, the actual errors are 


xa 


Y. 

ya, 


= Y 


Va' 


'cos ^ 


sin ip 


+ Va - g cos <p + a, sin <)>) 
S 


^-sin \p s 


cos \p 


+ HOT 1 s (3.16) 


The second term is lateral to the flightpath and consists of low-frequency transients 
when the DG or aircraft sideslip dynamics are excited; it is a bias in the case that 
a steady sideslip is imposed by the pilot. The filter interprets these errors as an 
equivalent crosswind in the amount 

W • j_ = Va(\j; - g cos <p + a sin cf> ) (5.8) 

L g 

Sideslip dynamics are omitted from the aircraft motion simulation so that the cross- 
wind accuracy in figure 5.1, rms(Wy), is somewhat optimistic. However, figure 5.5 
illustrates the effect of sideslip; the actual winds are null and sideslip is a 10° 
square wave. As seen, the ensemble mean crosswind estimate, E[Wy] , rises to the 
value of VaB predicted from equation (5.8) (9.6 m/sec), and the longitudinal wind 
estimate is unaffected. Additional effects include transients in lateral velocity 
and acceleration, y , y, excited by the brief, impulse-like excursions in B which 
the filter misinterprets as nonzero crosswind shear. 

Thus, the air velocity measurements detect a combination of crosswind and B> 
but they are individually unobservable. The addition of an air velocity vector 
sensor (e.g., refs. 31 and 32) or B-vane is necessary for their observability and 
the control of B or lateral air velocity. Aircraft have natural weather vane sta- 
bility so that control of B is not necessary for attitude stabilization. However, 
if an estimate of B is calculated from the filter output and fed back for attitude 
control, then the combined estimator and control system is neutrally stable in B* 

Such an estimate can be derived, using the theoretical relation 

sin B = (i^ip^ (5.9) 

where is the lateral body axis ((2b ) r is the second row of T^ r in table 2.2) 

and ip a is the air-velocity direction (along V - W) . Since zero sideslip was 
assumed in the filter, the estimate of B from equation (5.9) and the filter output 
must also give zero, approximately; that is, 

sin B = <1^ > i pa > 3 0 

and, therefore, 

\ 

3 = 3 

The transient response of the filter to sample wind histories is examined in 
figure 5.6. Figure 5.6(a) shows response to a mean wind sample with short periods of 
high shear (1 m/sec 2 ); response time for the average estimates, E[W x ],E[Wy] is simi- 
lar, using either VORTAC or MODILS data and on either axis. Excursions in rms errors 
associated with the shear events appear in the velocity and accelerometer bias states. 
These effects result from over-weighting acceleration bias and prior velocity estima- 
tion error as the causes of the higher-than-expected airspeed residuals which occur 
during periods of higher-than-expected wind rates. As seen in figure 5.6(a), the 
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effect on the longitudinal component of a., 
[rms(x)], is pronounced whereas the lateral com- 
ponent shows little sensitivity. These effects 
are also present during periods of higher wind 
rates in a turbulence sample. This observed 
influence of wind variation on estimation errors 
for the inertial states, V r , a r , is suppressed 
if an inertial grade IMU is used, since the 
weighting of acceleration and prior inertial 
velocity estimation error as the causes of air- 
speed residuals is much smaller relative to 
the prior wind-estimation error. 

Figure 5.6(b) shows response to a single 
turbulence sample with higher-than-average 
intensity. Here, wind shear is nonzero almost 
everywhere. Root-mean-square wind estimation 
errors resemble those of figure 5.6(a), with 
excursions superposed in association with wind 
shear. These errors are about the same size as 
those of the turbulence sample itself. 

Although rms velocity and acceleration errors 
are excited by wind variations, they are seen to 
be much smoother than the rms wind errors and 
are not significantly degraded from the results 
in figure 5.6(a); that is, accuracy for these 
states is only moderately affected by the tur- 
bulent wind variations. 





The filter models the wind with linear, 
first-order state equations forced by white 
noise and the parameters (i w ,a w ) whose values 
are to be selected. This only approximates the 
actual wind, which is simulated as a superposi- 
tion of independent constant and turbulence com- 
ponents. The filter’s parameters (x w ,a w ) have 
been selected at 20 sec, corresponding to the 
more rapid average rate of the turbulence, and 
at 6 m/sec, corresponding to the larger 
standard deviation of the mean wind. This 
choice results in the insensitivity of the 
velocity and acceleration bias estimate to the 
turbulent wind variations noted above and the 
relatively good agreement of the mean estimated 
wind with the total wind. 


Dead Reckoning 

In practice, position data are occasionally 
lost for periods of varying length, during which 
time position estimation errors diverge. For 
unaided dead reckoning, error growth is given by 



Figure 5.5.- Effect of steady side- 
slip on wind estimation. 
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(a) RESPONSE TO MEAN WIND VARIATION 


(b) RESPONSE TO TURBULENCE 



In the case in which am(t) is constant, velocity and position errors diverge linearly 
and quadratically with time, respectively. The amount of divergence for a given data 
dropout interval depends on the initial errors and on the IMU accuracy. The addition 
of airspeed aiding stabilizes the velocity error, and position error diverges only 
linearly. 
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The potential magnitude of these effects 
during 100 sec of dead reckoning in a maneuv- 
ering segment of the STOL approach outside 
MODILS coverage is shown in figure 5.7. 

These are unfavorable conditions for initial 
errors and IMU accuracy, but ones that none- 
theless occur in practice. Without airspeed 
aiding, velocity errors diverge to 15 to 
25 m/sec and position errors reach 1,000 m. 
With airspeed aiding, accuracy is clearly 
improved; velocity errors are stabilized, but 
position errors diverge by 500 to 600 m after 
100 sec, and significant maneuvering to 
return to the reference trajectory is still 
required when the reception of position data 
is resumed. Naturally, a high accuracy IMU 
diverges much more slowly, as shown in fig- 
ure 5.7, and provides for much longer 
periods of dead reckoning before errors 
become large. Here, the divergence is negli- 
gible after 100 sec. 


Vertical Axis Filter States 

An important design issue in this study 
was the trade-off between computation time 
and accuracy associated with the subdivision 
of the 14-state system into an independent 
4-state subsystem associated with vertical 
motion, (b^,z,z,z), and the remaining 
10 states associated with motion in the 
horizontal plane. The computation time for 
multiplying system matrices rises with the 
cube of system order; the proposed partition- 
ing would, therefore, reduce this time by 60% 
and result in a corresponding major reduction 
in the time required to execute the filter. 

The dynamics of the two subsystems 
(table 3.2) are mutually independent. In 
addition, the navaids can be similarly 
separated approximately. After partitioning 
the state into lower-order vectors associated 
with the two subsystems, denoted by 



then the dependence of the measurement 
residuals on these states can also be 
partitioned. 


STANDARD FILTER 

WIND STATES DELETED 

HIGH ACCURACY IMU 

(a = 0.001 g) 



POSITION 

NAVAID 


VORTAC 


DEAD RECKONING 


Figure 5.7.- Performance during dead 
reckoning. 
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y = Vh + Vv + * 


and the separation of the navaids can be made, provided the residuals for every 
navaid depend statistically almost entirely on the error state of one or the other 
subsystem; that is, if the residuals satisfy 

E[ (h v X y ) 2 ] « ElOi^n) 2 ] 
or 

E t (hyX v ) 2 ] » E[( W 21 

The gradient vectors for the navaids have been partitioned in this manner in 
table 3.4, and a review of that table shows that they divide into two sets: resid- 
uals for the derived MODILS altitude and altimeters are independent of and 

residuals for VORTAC, MODILS azimuth and. range, and air velocity have only second- 
order or smaller gradients with respect to the vertical-axis states, compared to 
their sensitivities to Xjj. This assumes that the elevation angles, EL m ,EL t are 
small: EL m is always small because of its coverage restrictions (table 2.3) and 

ELfc can reach high values temporarily on some trajectories. However, the dependence 
of VORTAC residuals on the vertical position error remains statistically negligible 
since o z << a x ,a y . 

This separation of navaids degrades the accuracy of both subsystems by neglect- 
ing some information on one subsystem while erroneously assigning the entire residual 
to the error states of the other. The actual accuracy loss was evaluated in simula- 
tion tests, using the reference STOL approach trajectory, and was found to be negli- 
gible for all states; therefore, the independent filter formulation is adopted in 
this study. 

Accuracy for the states of the vertical axis filter on the straight line test 
path (fig. 5.1) shows several features. First, baroaltimeter bias accuracy is 
invariant. This occurs because this bias is unobservable in the absence of altitude 
measurements from an independent source and will therefore be fixed on any approach 
until MODILS elevation data are acquired. Second, vertical position error is domi- 
nated by the baroaltimeter bias, and the accuracy is the same as for the bias; on 
sample cases with null bias, position accuracy is the same size as the altimeter 
noise (1.5 m) . Third, vertical velocity accuracy settles rapidly to 0.3 m/sec and 
does so independently of the random constants of each flight. This is considerably 
better accuracy than for the horizontal-plane velocity coordinates and results from 
the following factors: (1) the altitude-measurement bias error is constant and 

nearly uncorrelated with z; (2) the position measurement noise is much smaller for 
the altimeter than for VORTAC; and (3) the vertical acceleration measurement bias 
is nearly constant and can be accurately estimated from altimeter measurements. 

Note also that accuracy for all vertical axis states is independent of the position 
navaid used by the horizontal plane filter. Further, accuracy is insensitive 
to vertical-axis maneuvering (flares and pitchovers) because of (1) accurate calibra- 
tion of the vertical acceleration measurement and (2) the low values of |z|, below 
0.10 g, used in passenger operations. These favorable factors are present on all 
flightpaths appropriate for passenger operations, and the good, invariant accuracy 
seen here will be achieved throughout the terminal area. The vertical wind, W z , is 
not a filter state; it is unobservable but consists of turbulence only and does not 
affect accuracy for the vertical-axis inertial motion. 


62 



Flightpath Angle and Angle-of-Attack Estimation 

Inertial flightpath angle and angle of attack are variables of interest in auto- 
matic control systems and can be derived from the estimator's output with an accuracy 
that depends largely on the performance of the vertical axis filter. Formulas for 
estimating these angles are 


y = sin 1 (-i/V) 
a = sin -1 ( (k^ ».ip a > / cos (3) 

and expressions for estimation errors, derived from equation (5.10) are 

Y = <k r ,V>/V + HOT's 

a ■ %-ira > + < ib'ip a > 

where 


-Pa 


-k y + i n ik T 
—Pa a -Pa Va 


Y, = 


-<k Pa .V - W>/Va 


^Va = <1 p ~ W>/Va 


>(5.11) 


MPa 


'sin <f> cos ip - cos <J> sin 0 sin \p 

(kj J ) r = - ( J lj 3 ) r $ + (i^) r cos <|>0 + I sin <f> sin iJj - cos <f> sin 0 cos ] ip 


g 


Expressions for the runway axes components of the vectors ib > k^ »ip a » J-Pa ’iiPa ’ 
required to evaluate the inner products above are given by the rows of Ibr’^p r 
(table~2.2) and can be computed from <|),0,iJj,V r ,W r . The errors y>« depend on 3 
V r ,W r ,* g ,0 variational expressions derived from equation (5.11), assuming that 

YjQj'I'Va “ v are small angles, are 


Y = - # + HOT's 

a = cos (f>[0 + (I - W ) /Va] + HOT's 

8 z j 


(5.12) 


As seen, y depends on z and inversely on inertial speed; a has several significant 
error sources (§g,£,W z ) and varies inversely with airspeed. 


Ensemble averages for y and a on a straight line path are shown in figure 5.8; 
speed is varied, over the approach speed range (60 to 140 knots) during the test, as 
plotted in the figure. Root-mean-square (z) is invariant at 0.3 m/sec and rms(Y) 
varies inversely with speed, as expected, from 0.25° to 0.6°. Results for rms(a) 
show the expected dependencies on 0 g ,W z ,I, and Va; accuracy is much poorer than for 
Y and varies from 0.7° in favorable conditions (higher speed, no pitch error) to 3° 
in unfavorable conditions (low speed, large pitch error). The dominant error sources 
are the unknown turbulence, W z , with a time average of 0.7 m/sec here, and the VG 
pitch error. These results imply significant corresponding errors in estimating or 
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controlling lift force directly; for conven- 

tional aircraft (l/W)(dL/da) = Cl^/Cl, and this 

varied from 0.04 to 0.2 g/deg over the approach 
speed range for an example aircraft. 

Vertical Axis Filter Parameter Values 

The performance effects of varying the 
selected values of the filter’s state process 
noise correlation time can be examined readily 
for the vertical axis filter. These parameters, 
T 2 , T b^, represent the time scale at which 

acceleration measurement and altimeter biases 
are varying. They control the rate at which 
the accuracies of the estimated biases are 
assumed by the filter to degrade in the absence 
of new measurements and the weighting given to 
information from new measurements relative to 
the prior information. In the present simula- 
tion, the two measurement biases are constant, 
or nearly so, and the filter is easily matched 
to the simulation model by selecting the same 
correlation times as in table 3.2 
( T z ,T hb ~ 10 3 »10 4 sec). The effect of various 
mismatches of the filter’s values from those 
values are shown in figure 5.9; results for 
four cases with the following parameter values 
are given: 


Case 

1 

2 

3 

4 


10 : 

20 

10 3 

20 


T hb 

10 4 

10 4 

50 

50 


TIME, sec 


Figure 5.8.- Accuracy of inertial 
flightpath angle and angle of 
attack. 


The figure includes a sample case history along 
with the theoretical standard deviation from the 
filter’s covariance and the ensemble average. 

The existing design values (case 1) result in 
good agreement between theoretical and actual 
accuracy in both z and z and show the best 
transient and steady state performance of the 
four cases. If x~ is reduced (case 2), excess 
weight is placed on acceleration biases as the 
cause of measurement residuals, and the results show the poorest performance among the 
four cases with degraded theoretical and actual steady state performance, particularly 
in acceleration. A reduction of T^b (case 3) reduces the capacity of the altimeter 
measurements to calibrate the acceleration bias and the results show degraded theoreti- 
cal transient response but retain good actual steady state accuracy. Finally, case 4 
has reduced values for both correlation times; the filter predicts poor accuracy for 
both z and z, with little capacity to calibrate the acceleration measurement or 
improve the initial accuracy (if the selected correlation times were accurate it would 
be pointless to process the measurements in this case). However, the actual accuracy 
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is similar to that of case 3; transient response is slower than case 1 but the 
acceleration bias is successfully estimated and steady state accuracy is good. 

These data illustrate some significant performance changes with the order of 
magnitude of errors in the selected values of T z ,T bj 1 ^* Inaccuracies affect perfor- 
mance in one or more states and can result in reduced transient response and steady 
state accuracy and in increased disagreement between predicted and actual accuracies. 
Therefore, if the actual measurement bias processes differ significantly from the 
simulated ones, parameter adjustments based on recorded flight data would be useful in 
optimizing filter accuracy. 


Summary 

Performance details have been studied using simulation results from a straight 
path near the runway, and other data. Several general conclusions on filter design 
issues and estimation accuracy trends are indicated. 

1. The 14 states can be separated into independent lower-order systems asso- 
ciated with vertical axis motion (4 states) and motion in the horizontal plane 

(10 states) with negligible loss of estimation accuracy for any state and a large 
reduction in the required computation time to process a scalar measurement. 

2. The wind states and air velocity measurements cannot be separated from the 
horizontal plane filter into an independent wind filter without a significant loss of 
accuracy in position, velocity, and wind estimates during VORTAC use and dead 
reckoning. 

3. Selection of parameter values for the state noise and measurement error 
models is based on simulation tests with realistic models or recorded flight data and 
is generally aimed at a reasonable fit of the observed state and measurement error 
variations, at reasonable agreement between Monte Carlo rms errors and the filter* s 
covariance for those states sensitive to the parameter values, and at minimizing the 
resulting ensemble averages for these states. 

4. Horizontal plane position estimation accuracy depends principally on position 
measurement accuracy; IMU accuracy improvements yield only limited position-accuracy 
improvements. Before entering MODILS coverage, position accuracy is dominated by 
VORTAC calibration accuracy. Although information on the VORTAC biases is theoreti- 
cally available from the IMU and by cross-calibration, it was found that this infor- 
mation was negligible until close to the runway, when passing near the VORTAC station 
permits range bias calibration to about 60 m by the IMU, and entry to MODILS coverage 
permits calibration of the VORTAC biases to the accuracy of the MODILS measurements 
(~15 m) . An inertial grade IMU would increase range calibration accuracy when pass- 
ing close to the VORTAC station and would provide a capacity to calibrate bearing to 
better than 0.5°, if a turn were executed for this purpose. 

5. The accuracy of horizontal plane velocity estimation is sensitive to both 
position measurement and IMU accuracy; improvements in both types of sensors are 
required over the present sensors to achieve 1 m/sec accuracy throughout the terminal 
area. During VORTAC use, lateral velocity is biased proportionally to bearing bias, 
and noticeable velocity and wind-error transients occur when passing near the VORTAC 
in samples with larger-than-expected range bias. 
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6* The accuracy of horizontal plane wind estimation is dominated by inertial 
velocity estimation accuracy throughout VORTAC use, and rms wind estimation errors 
are about the same size as the turbulent component of the wind. In addition, cross- 
wind is indistinguishable from sideslip and directional gyroscope error so that cross- 
wind and 3 are unobservable and 3 is uncontrollable. Vertical wind and a are 
also unobservable, and additional sensors are necessary for the observability of 
a and 3> and of the direction of .he air velocity vector. 

7. Vertical axis inertial motion (z,2,z) is generally estimated to much greater 
accuracy than the horizontal plane motion (30 m, 0.3 m/sec, less than 0.01 g) . This 
accuracy or better is achieved throughout the terminal area. 

8. In dead reckoning, airspeed measurements stabilize the velocity estimation 
error, but rms position error rates of 5 m/sec can occur on maneuvering flightpath 
segments, and significant maneuvering to return to the reference path is required 
after 100 sec of dead reckoning. An inertial grade IMU would provide at least an 
order of magnitude improvement. 


6. ESTIMATION ACCURACY ON A V/STOL APPROACH 


Accuracy during the reference STOL approach trajectory (fig. 2.1) is evaluated 
in this section. The filter is initialized at approximately the terminal area entry 
distance (25 n. mi.). Results are shown for the initial 200 sec, to illustrate 
initial transients, and for the last 500 sec, to illustrate accuracy near the runway 
and during maneuvering. The results obtained are typical of most approach paths. 
Accuracy is shown graphically (fig. 6.1) and as a table of values for various points 
along the approach (table 6.1), and its general behavior for all states of interest 
in the control is summarized in table 6.2. In addition, a sample case and the ensem- 
ble extremes are reviewed to indicate the range of error histories encountered in an 
approach (fig. 6.2). Scales are changed in these figures where useful for greater 
resolution. Accuracy is given for the path axis components of the kinematic states; 
these are useful working axes for an automatic control system, which is assumed 
decoupled in these axes so that system exitation owing to estimation errors is 
similarly decoupled. 


Measurement Biases and Navaid Selection 

Results for the VORTAC biases confirm the conclusions of the previous section: 
these biases are poorly observable until passing close to the VORTAC station 
(500 sec < t < 550 sec) at high bearing rates (3° /sec), and favorable IMU accuracy 
(0.01 g) results in range bias calibration to 30 m by the IMU. Bearing bias is very 
poorly observable until entry to MODILS coverage (at t = 740 sec) permits VORTAC 
calibration to about 20 m and 0.2°. During MODILS coverage, the VORTAC is well cali- 
brated, but its equivalent position fix accuracy remains much poorer than that of 
MODILS because of signal noise; therefore, VORTAC data processing can be dropped with 
negligible loss of accuracy. VORTAC processing resumes in the event of a failure 
of MODILS reception or of exit from MODILS coverage in a missed approach. At these 
events, both of which occur in practice, the VORTAC calibration is renewed and its 
good accuracy avoids the maneuvering to a biased position estimate that would arise 
from large calibration errors. 
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POSITION ACCURACY 









(b) Velocity and winds. 


Figure 6.1.- Continued 
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Figure 6.1.- Concluded. 
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TABLE 6.1- PERFORMANCE SUMMARY — ENSEMBLE RMS TRANSLATIONAL-STATE ESTIMATION ERRORS ON A STOL APPROACH 



Location on 
reference 
trajectory, m 

Navaid a 

Position 
errors, rms 
m 

9 

Velocity 
errors , rms 
m/sec 


Acceleration 
errors , rms , 
g 


X 

y 

z 


Longi- 

tudinal 

Lateral 

Normal 

Longi- 

tudinal 

Lateral 

Normal 

Longi- 

tudinal 

Lateral 

Normal 

Initialization 

629 

45730 

-1161 

V.B 

245 

1475 

27 


7.8 

0 

0.012 

1 


First segment 

629 

31300 

-1161 

V.B 

256 

1260 

29 


2.6 

0.3 



■ 

First turn entry 

629 

-915 

-1161 

V.B 

47 

55 

29 

1.8 

3.5 

.3 

.01 

.007 

.007 

First turn exit 

-895 

-2438 

-1036 

V.B 

75 

30 

30 

1.4 

4.0 

.3 


wwm 

.006 

Second turn entry 

-4298 

-2438 

-1036 

V.B 

64 

115 

30 

1.4 

2.0 

.3 


iBH 

.006 

Second turn exti 

-4298 

0 

-1036 

M,B,E 

10 

14 

5 

1.2 

1.2 

.3 

.013 

.006 

.005 

Midhelix 
Glide slope 

-2807 

-2440 

-627 

M,B 

14 

12 

4 

2.0 

1.2 

.5 

.025 

.025 

.005 

capture 

-2807 

0 

-318 

M.B.E 

11 

12 

3 

1.8 

1.0 

.3 

.020 

.012 

.005 

Glide slope 

-1500 

0 

-207 

M,B,E 

11 

8 

2.5 

1.2 

.6 

.3 

.010 

.007 

.005 

Glide slope 

-750 

0 

-109 

M,B,E 

10 

6 

2.5 

1.2 

.6 

.3 

.010 

.007 

.005 

Flare start 

-236 

0 

36 

M,R 

14 

4,5 

2.5 

1.0 

.6 

.3 

.010 

.007 

.005 

Landing zone 

>-100 

0 

<10 

M,R 

10 

3 

1 

1.0 

.6 

.3 

.010 

.007 

.005 


V - VORTAC, M - Modils Az and DME, B = Baroaltimeter, E = Modils elevation, R = Radar altimeter. 























TABLE 6.2.- SUMMARY OF ESTIMATION ERROR CHARACTERISTICS FOR A TERMINAL-AREA 

APPROACH TRAJECTORY 


State 

Error behavior 

Significant 

error 

sources 

Range of 
rms errors 

Position 

Longitudinal 

Heading, location dependent; typically 
piecewise stationary 

VORTAC range 
MODILS range 

10-250 m 

Lateral 

Heading, location dependent; rms 
typically linear in range 

VORTAC bear- 
ing, MODILS 
azimuth 

6-1500 m 

Normal 

Invariant or linear rms 

Altimeters, 

elevation 

1-30 m 

Velocity 

Longitudinal 

Maneuver dependent; stationary in 
equilibrium flight 

Heading, 

pitch, 

accelerometers 

1.5 to 
5 m/sec 

Lateral 

Maneuver dependent; biased by VORTAC 
bearing bias 

Roll angle, 
VORTAC biases 

0.6 to 
5 m/sec 

Normal 

Stationary 


0.3 m/sec 

Acceleration 

Longitudinal 

Maneuver dependent 

Heading, 

pitch, 

accelerometers 

0.01-0.04 g 

Lateral 

Maneuver dependent 

Roll, 

accelerometers 

0.01-0.03 g 

Normal 

Stationary 


0.005 g 

Attitude 

Pitch, roll 

Maneuver dependent 

VG 

to 2° 

Heading 

Maneuver dependent 

DG 

to 6° 

Wind 

Longitudinal 

Contains V 

V 

1-4 m/sec 

Lateral 

Contains V, indistinguishable from 
3, ^g. Slow transient response. 

V,g,DG 

1-6 m/sec 

Vertical 

Not estimated 

Unobservable 


System in dead 
reckoning 

Linear position divergence 

IMU 

Drift to 
5 m/sec 
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(a) Measurement biases and positions. 


Figure 6.2.- 


Sample case and extreme estimation errors. 
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The baroaltimeter bias is unobservable in the absence of an independent alti- 
tude measurement, so that the bias error is stationary at 30-m accuracy until the 
nearly unbiased elevation data become available in midturn (at t = 785 sec) and 
improves calibration accuracy an order of magnitude to 3 m. Subsequently, elevation 
data drops out temporarily in the helix (825 sec < t < 925 sec) because of heading 
limits on signal reception, but reception resumes near the start of the glide 
slope, after which the derived altitude measurement accuracy of the elevation 
data and the corresponding baroaltimeter calibration improves and reaches 1 m at 
landing . 


Baroaltimeter data are processed throughout the approach up to the landing zone, 
whether or not elevation data are available because the baroaltimeter measurement has 
lower noise standard deviation (o v = 1.5 m) than the derived elevation altitude 
(6.5 m at entry to elevation coverage on the present test path) and yields a more 
accurate vertical velocity estimate than do elevation data alone. This is readily 
demonstrated in simulation tests (not shown). Conversely, elevation data are nearly 
unbiased and yield a more accurate altitude estimate throughout the coverage of alti- 
tude than does the baroaltimeter; results for the remaining states of the vertical 
axis filter show that this is the principal benefit of processing elevation data. 

The baroaltimeter bias is constant in the present simulation model so that no 
loss of calibration accuracy occurs in figure 6.1, once the calibrating signal is 
lost. However, this result is optimistic; flight data show a significant variation 
in baroaltimeter bias with altitude over the altitude range of this test so that the 
altimeter calibration cannot be maintained in practice, following the loss of the 
calibrating signal. Additional effects are that (1) altitude estimation accuracy will 
degrade more rapidly toward its a priori value in the absence of new calibrating data, 
(2) vertical velocity will degrade slightly on nonlevel paths because the fictitious 
vert i- ca l velocity, h(dbhb/dh), implied by altimeter bias variations or because of 
calibration lags on the glide slope, and (3) the appropriate time constant for the 
filter's baroaltimeter bias process noise model is smaller than that selected here. 
Some possibilities for improving this situation are (1) redesign of the receiving sys- 
tem to remove the heading restriction on elevation signal reception, and (2) a study 
of actual baroaltimeter bias errors to determine if a more accurate model would permit 
on-line determination of their altitude dependence to a useful degree. 


A 

b 


hb 


Finally, we note that it is possible 
to simplify the vertical filter by remov- 
ing the baroaltimeter bias state to a 
separate filter where it is calibrated by 
the elevation data (see sketch H) of an 
independent baroaltimeter calibration fil- 
ter. However, this trades a small reduc- 
tion in the computation time of the com- 
plete filter algorithm (about 3%) for some 
loss in altitude accuracy during elevation 
coverage. Similar separate filters to 

calibrate VORTAC biases from MODILS data are also possible, but they result in sig— 
nificant accuracy losses outside MODILS coverage. These devices could, however, be 
useful in applications in which lower accuracy is acceptable and as a means of adding 
the capability for on-line calibration to the conventional complementary filter 
estimator. 
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Position 


Longitudinal accuracy is dominated by range-measurement accuracy for most of the 
approach, which is nearly along a VORTAC radial until close to the runway, and along 
a MODILS radial during most of the final segments. This approximate relation is 
expected to be typical of approach paths for the present antenna sitings. Thus, 
accuracy (fig. 6.1 and table 6.1) settles to a nearly constant value at 225 m until 
near the VORTAC station, reaches 30 m at closest approach to the station (t = 640 sec), 
and further improves to 10 m after entry to MODILS coverage; it then remains there 
until landing. This accuracy sequence reflects the VORTAC range calibration accuracy 
history and the invariant MODILS range measurement accuracy. Briefly, longitudinal 
errors are approximately piecewise stationary and change by more than an order of 
magnitude during the approach. 

Lateral position accuracy depends primarily on the VORTAC bearing and MODILS 
azimuth accuracies, which results in range dependent lateral position accuracy. Ini- 
tial accuracy is poor at 1,500 m and drops linearly with range to the MODILS station 
during the glide slope segment, reaching 3 m at landing. Lateral position accuracy 
changes by three orders of magnitude during an approach, with periods of approximately 
linear, step, or invariant behavior. 

Normal position accuracy follows the baroaltimeter accuracy closely; it is 30 m 
until elevation data are available, after which it improves rapidly to 5 m and then 
declines on the glide slope, reaching 1 m at landing. The accuracy change exceeds 
one order of magnitude for the approach, with periods of invariant, step, or linear 
behavior. Finally, it is apparent in figure 6.1 that position accuracy varies widely 
between the three control axes (errors are nonisotropic) , with order-of-magnitude 
differences in some segments, particularly at entry to the terminal area. 


Velocity, Wind, and Acceleration 

Initially, longitudinal velocity accuracy settles to a steady state value 
(1.5 m/sec) in about 30 sec; this repeats the initial accuracy behavior previously 
seen in the straight-test-path results, because the relevant measurements (range, 
acceleration, airspeed) and their accuracies are identical. Subsequently, accuracy 
shows modest excursions to 2 to 4 m/sec in the half turn (legs 4 and 5) and helix 
(leg 7) and settles to 1 m/sec for the glide slope and landing zone. The excursions 
can be traced to excursions in \pg and corresponding longitudinal acceleration 
measurement error transients that recur in most samples. 

Initially, lateral velocity accuracy settles to 3 m/sec; the settling transient 
is much slower and the steady state accuracy poorer than for the longitudinal axis as 
a result of (1) the larger position measurement noise and lack of lateral air velocity 
information, and (2) the presence of a bias error proportional to the bearing cali- 
bration error. During MODILS coverage accuracy improves gradually to 1 m/sec and 
then shows excursions, because of maneuvering, followed by progressive improvement 
on the glide slope with improving position and acceleration measurement accuracies, 
reaching 0.6 m/sec at landing. 

Normal axis velocity accuracy settles rapidly to 0.25 to 0.5 m/sec, and this is 
obtained uniformly throughout the approach independent of location or the availability 
of elevation data. 
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Longitudinal wind accuracy is dominated by the inertial velocity error and 
se itles rapidly to 1 to 2 m/sec for most of the approach, with excursions during 
legs 4, 5, and 7 paralleling those of the inertial velocity. Crosswind accuracy 
follows the lateral inertial velocity accuracy closely, with excursions resulting 
from ip g excursions superposed (legs 4 and 5, particularly). Thus, the DG error 
transients during turns appear in both the longitudinal and lateral wind estimates 
via different routes. Further, the crosswind results are modestly optimistic, since 
sideslip also appears in its errors, and sideslip dynamics are omitted from the motion 
simulation. 

The acceleration accuracy history is consistent with the conclusions of the pre- 
vious sections: accuracy exceeds 0.01 g on the unaccelerated segments (throughout 

the initial straight leg and the glide slope) and for the normal axis over the entire 
approach, but degrades in maneuvering segments to 0.025 g to 0.04 g for the longi- 
tudinal and lateral axes. 


Accuracy Summary 


The estimation accuracy encountered in this example STOL approach is expected to 
be typical of most terminal area approach trajectories with similar navaid sitings, 
both in magnitude and behavior. The empirical results for various points along the 
approach and the observed trends in error behavior and the principal error sources 
for all states used in the control are recapitulated in tables 6.1 and 6.2. 

Table 6.2 indicates that the system of this study is characterized by complex naviga- 
tion errors along a terminal area approach. Accuracy is best for the normal axis 
motion, with simple behavior for position and invariant velocity and acceleration 
accuracy independent of location, maneuvering, and navaid. For the horizontal plane 
motion, position accuracy tends to be piecewise linear or stationary, with order of 
magnitude changes along the approach and differences among control axes, and velocity 
and acceleration accuracy excursions occur during speed changes and turns. 


Results with an Inertial Grade IMU 


Figure 6.1 includes a comparison of results for the present system with accura- 
cies achieved using an inertial grade IMU (0.001 g, 0.01° accuracy for acceleration 
and attitude measurements, respectively) in order to show the performance improve- 
ments and its limitations from this sensor improvement. The results show a somewhat 
sstllsr and more accurate calibration of the VORTAC range when passing near the 
VORTAC station and exhibits the capacity of a sufficiently accurate IMU to calibrate 
bearing (to 0.2 ) during the quarter turn (leg 2). The position accuracy shows no 
important gains, but we note from earlier analyses that (1) a turn executed early in 
the approach would provide a significant improvement in lateral accuracy throughout 
VORTAC use from the IMU's calibration of bearing, and (2) position accuracy in dead 
reckoning would improve to an important degree. 

The velocity components show major improvements in accuracy. The velocity error 
excursions during maneuvering owing to the pendulous attitude gyrodynamics of the 
standard IMU have been eliminated, and accuracy is well below 1 m/sec for the longi- 
tudinal axis over nearly the entire approach and for the lateral axis during MODILS 
coverage. On the lateral axis, the slow initial settling time cannot be reduced 
significantly by the improved IMU, nor is the bias caused by bearing calibration 
error significantly reduced (although it can be by executing a turn early in the 
^PP^o^^-L) . On the longitudinal axis, the modest excursion in accuracy during the 
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quarter turn is due to mapping of the remaining lateral velocity bias onto the longi- 
tudinal axis as the path axes rotate. Longitudinal wind accuracy is improved to the 
accuracy of the airspeed measurement nearly everywhere. Crosswind accuracy is also 
improved by the absence of DG error transients and throughout MODILS coverage; the 
effect of lateral velocity bias during VORTAC use is unchanged here (but can be nearly 
eliminated), but the lack of observability of crosswind, sideslip, and vertical wind 
is unaffected by IMU accuracy. 


Sample Case and Ensemble Extreme Errors 

Error histories for a sample approach, together with an envelope of extreme 
errors for the ensemble of 10 approaches, are shown in figure 6.2. Estimation errors 
are not ergodic; that is, sample case time averages can differ considerably from the 
ensemble averages. Figure 6.2 provides a view of typical error histories experienced 
by the control system during an approach and the range of these errors. The ensemble 
extreme magnitudes are generally twice the rms values for all states at all times; 
that is, 2a errors are experienced within 10 approaches in which all random varia- 
bles are randomly sampled. 

The VORTAC range bias error sample follows the behavior of the ensemble average 
closely; that is, the error is nearly unobservable and constant until passing near 
the VORTAC station results in an error reduction to 30 m and until a further reduction 
to 15 m occurs at entry to MODILS coverage. Other samples with significant range bias 
show the same calibration error behavior. The VORTAC bearing bias error sample shows 
an initial loss of accuracy by 1° as the filter settles the combined initial lateral, 
velocity error and bearing bias error to a combination that satisfies equation (5.6); 
subsequently, the error remains constant until MODILS calibrates bearing to 0.1 . 

The same settling to a combination that satisfies equation (5.6) occurs on all samples 
after initialization, but the steady state error can be fortuitously reduced as well 
as increased; one sample showed a reduction from 2° bias error to nearly perfect cali- 
bration. The ensemble of 10 bearing bias samples (table 4.1) is significantly skewed 
from the statistical distribution of b tb and contains only one sample with 
b tb > 0; this sample is the one mentioned immediately above, which settles to zero 
steady state error, and which defines the unexpected ensemble maximum and nonzero 
mean extreme seen in figure 6.2(a). Nevertheless, bearing bias is nearly unobserv- 
able in this segment of the approach, and a larger ensemble would give the expected 
invariant extremes and zero ensemble mean value. Corresponding effects appear in the 
extremes of the lateral position and velocity errors discussed below, and similar 
effects appear in the baroaltimeter bias and normal position error extremes as a 
result of skewing of the b bb ensemble. The baroaltimeter bias error is unobservable 
and constant in the absence of elevation data, as seen in both the sample case and 
extremes. Otherwise, the extreme magnitude is reduced to 6 m after entry to elevation 
coverage; this improves to 2.5 m on the glide slope and to 1.5 m at landing, where 
the radar altimeter is used. 

The longitudinal position error sample is biased by the VORTAC range bias 
initially with smaller variations superposed because of VORTAC noise (these errors 
are of the order of 50 to 75 m in samples with no range bias). On the glide slope, 
the sample shows a mean negative value, which results from an unusually large sample 
value of the MODILS range bias (-12.8 m) . Lateral position error sample is initially 
biased by the VORTAC bearing bias errors; it degrades initially with b tb and then 
improves linearly with reducing range to the VORTAC. On the glide slope, it is 
unbiased because of the small sample case MODILS azimuth bias. However, bias is sta- 
tistically the dominant azimuth error, and in many other samples lateral position error 
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is dominated by a range-dependent bias in proportion to the azimuth signal bias. The 
normal position error sample is initially dominated by the baroaltimeter bias, with 
small variations superposed (±5 m) because of altimeter noise. After entering eleva- 
tion coverage, this error is within ±10 m, a result equally of the altimeter noise 
and of the remaining bias error, and then improves to ±5 m on the glide slope, where 
the bias error is much smaller. During the final approach, the ensemble extremes 
have properties similar to the ensemble averages; that is, they are relatively con- 
stant for longitudinal position (±20 m) , decrease with distance to the runway for 
lateral position (±20 m declining to ±5 m) , and are ±5 m for normal position, improv- 
ing to ±1.5 m at landing. 

The longitudinal velocity error sample is unbiased and within ±2.5 m/sec on the 
unaccelerated segments; during maneuvering segments, error excursions are induced by 
acceleration measurement error excursions, and the extreme magnitude expands to 
5 m/sec at several times. The lateral velocity sample is biased by the bearing bias 
error sample during VORTAC coverage but is unbiased on the glide slope, with extremes 
declining to ±1 m/sec at landing. The extreme magnitude is about 5 m/sec for most of 
the approach before the glide slope as a result of the effects of both VORTAC bearing 
bias errors and maneuver-induced lateral acceleration transients. The normal axis 
velocity error is stationary and unbiased and remains within 1 m/sec for the entire 
approach in this and all samples. 

As previously noted, acceleration errors are stationary and unbiased for the 
initial straight segment, on the glide slope, and for the normal axis over the 
entire approach; in these cases, the sample errors are typical of all samples, and the 
extreme errors are nearly constant at ±0.01 to ±0.015 g for all axes. During the 
maneuvering segments, the sample case longitudinal and lateral errors show increased 
low-frequency excursions, which are reflected in similar excursions or expansions of 
the ensemble extremes and an increase in the magnitude of the maximum error to 0.05 g 
at several times. 


7. DISCUSSION 


This discussion focuses on an evaluation of the adequacy of the present estima- 
tion system’s accuracy in applications to automatic and IFR operation of an aircraft 
in the terminal area along an area navigation route to a final approach and landing. 
One approach to the evaluation is a simulation demonstration that the estimation and 
control system satisfies the performance criteria on tracking errors and control 
activity that suffice for system acceptability. However, performance criteria for 
the flight operations of interest are incompletely known. Some safety criteria are 
given in FAA publications and other sources; they concern principally position track- 
ing accuracy, which is shown to depend mostly on low-frequency estimation errors. 
Criteria for ride quality and control activity, which depend principally on higher 
frequency estimation errors, are undefined. A complete evaluation is, therefore, 
beyond the scope of this work. Nevertheless, a generic automatic trajectory tracking 
system (4D) can be introduced into the simulation to determine the trajectory errors 
and control activity excited by the translational state estimation errors of the 
present system; moreover, the estimation accuracy required to satisfy the known 
operational criteria and the filter’s ability to meet these can be determined. 
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Effect of Estimation Errors on the Trajectory Tracking 
Performance of an Automatic Control System 

A simplified generic model of the combined automatic control logic and aircraft 
is derived in appendix D. This model does not require representation of any details 
of the aircraft force, moment, and thrust generation processes and suffices to deter- 
mine the effects of translational state estimation errors on control of the transla- 
tional degrees of freedom. Additional but independent effects arise from errors in 
the estimates of other states not considered here but required by the control system, 
such as those required to control attitude, engine, configuration, and actuator dynam- 
ics. Appendix D also contains analytical and simulation results, including transfer 
functions, frequency response amplitude ratios, and Monte Carlo evaluation of the 
output statistics along the STOL approach. 

Estimates of the translational states are used by the control to estimate the 
current trajectory tracking errors for feedback to the control laws, as outlined in 
appendix D. Estimation errors enter the control system through this feedback as 
differences between the true and estimated tracking errors; those differences excite 
erroneous control activity and corresponding tracking errors. The estimated path 
axes tracking errors can be written as 

6S p = T pr ( ^Va’V (zc r " * r > (7>1) 


where refers to any of the vectors R, V_, and ji, and where _zc_(t) is its current 

commanded value. The feedback control law is decoupled in path axes, and this results 
in the same decoupling of system response to estimation errors. The error in esti- 
mating the feedback quantities is 


6z = -T z + HOT* s 
p pr r 


( 7 . 2 ) 


Here, transformation angle errors introduce only second-order effects. Thus, the 
effects of estimation errors are defined by the control law acting directly on these 
errors and the rms estimation errors shown in figure 6.1 are also the rms input 
disturbances to the control. 


The analysis of frequency response amplitude ratios (appendix D) indicates sev- 
eral trends for the steady state response to sinusoidal estimation errors. Position 
dispersions respond principally to the low-frequency content of R, V, and a, that is, 
to frequencies below the control bandwidth (given as the natural frequency of the 
translational transient response dynamics imposed by the feedback control law). 

Biases in V or a do not appear in the velocity or acceleration dispersions but, 
instead, are converted by the control into position hang offs which are superposed 
with any bias in R. Velocity dispersions have peak sensitivity to the content of 
R, V, and a at the control bandwidth; that is, the control maps velocity estimation 
errors at this frequency into velocity dispersions of the same amplitude. Control 
activity can be calculated as the corrective acceleration required to impose the 
desired transient response on the translational states. For CTOL aircraft this 
acceleration maps into corresponding engine power, roll angle, and pitch angle activ- 
ity. It responds to R, V, and a principally at frequencies at and above the control 
bandwidth. The system response to estimation errors therefore depends on errors in 
all three states, R, V, and a and their distributions on the frequency spectrum 
relative to the control bandwidth, and on their mutual correlations. 


Simulation results for sample case and ensemble rms behavior with the present 
estimator (appendix D) indicate that position dispersions are predominantly at low 
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frequencies and that rms(6R) and rms(R) are the same size during most of the approach; 
however, there are significant excursions of rms(6R) in excess of rms(R) occurring 
during turns and in response to large changes in rms(R), such as the step changes 
seen at entry to MODILS and elevation coverage. Velocity dispersions are dominated 
by low-frequency transient variations caused by continual control effort to regulate 
position tracking errors, combined with the system response to V at frequencies 
near the control bandwidth. The result is that rms(6V) exceeds rms(V) by a factor 
of 1.5 to 2 in most circumstances. 

Control activity is principally at frequencies near and above the control band- 
width. Control authority for regulating tracking errors is typically 0.1 g in 
passenger operations, and it is desirable to limit the control activity owing to 
estimation errors in order to retain sufficient margin for regulating the effects of 
turbulence. For a given estimation-error environment, control activity increases 
with control bandwidth; the bandwidth can be selected for acceptable activity levels. 
The results, using suitable bandwidths, show rms activity that is roughly uniform 
over the approach for each axis; bandwidth and rms activity were 0.3 rps and 0.02 g, 
respectively, for the normal axis and 0.125 rps and 0.04 g for the longitudinal axis. 
In flight, turbulence maps principally into normal axis acceleration disturbances so 
that greater estimation accuracy is needed for this axis to achieve sufficient con- 
trol margin reserve and .bandwidth. For the longitudinal and lateral axes it is 
likely that the estimation errors of the present system are the dominant disturbances. 
For a given control bandwidth, control activity depends on sensor accuracies; a com- 
parison of activity for VORTAC and MODILS and for IMU f s of different accuracies indi- 
cates that control activity is dominated by position sensor accuracy when IMU accuracy 
is good, and is dominated by IMU accuracy when position accuracy is good. Thus, 
activity depends on both position sensor and IMU accuracies. It was found that an 
inertial grade IMU, combined with position accuracy equivalent to that of MODILS, was 
needed to obtain negligible control activity (0.01 g) because of estimation errors at 
the control bandwidths studied. 


Estimation Accuracy Requirements for Terminal Area Operations 

Some tracking accuracy criteria for various operations is available from FAA 
publications and other sources. Corresponding estimation accuracy requirements can 
be derived from these and compared with the accuracy achieved by the present estima- 
tion system. These tracking accuracy criteria generally reflect the position accuracy 
needed for a sufficiently low risk of collision or unsuccessful landing within the 
existing or prospective traffic separation standards and runway dimensions. In turn, 
these separation standards reflect the generally available or expected navigation aid 
accuracies relevant to the National Airspace System and its users (that is, accuracies 
corresponding to VOR/DME, ILS localizer or MLS, altimetry, and air data systems). 

Table 7.1 summarizes the maximum allowed standard deviations of the system dis- 
persions. Automatic landing requirements for CTOL indicate allowed dispersions of 
the touchdown point under Category II visibility conditions (ref. 33); those for the 
Space Shuttle Vehicle (SSV) are for dispersions at the nominal touchdown time on a 
45 x 3,050-m runway (ref. 34), and those for STOL operations are estimated assuming 
a 46 x 610-m runway. The requirements for CTOL Category II automatic approaches were 
derived from data in reference 35. In addition, criteria for IFR area navigation 
operations are defined in reference 36. These are compatible with use of existing 
VOR/DME equipment and conventional altimetry so that it can be assumed that the 
present system satisfies these operational criteria. 
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TABLE 7.1- SOME COMPARISONS OF ESTIMATION ACCURACY WITH OPERATIONAL REQUIREMENTS 


State 

Trajectory-tracking requirements, la 

Navigation accuracy requirements 

, lo 

Kalmal-filter 
estimation accuracy, la 

MLS 

RAINPAL 

Automatic landing 

CTOL, Cat. II 

Automatic 

landing 

CTOL, Cat. II 

Landing 

zone 

CTOL, Cat. II 

Cat. IV 
specifi- 
cation 

Landing 

zone 

CTOL 

ssv 

STOL 

Decision 

height 

At 

5 n . mi . 

CTOL 

SSV 

STOL 

VTOL 

Decision 

height 

At 

5 n. mi. 

Decision 

height 

At 

5 n. mi. 

Position, m 

















Longitudinal 

229 

73 

31 

* 

* 

132 

42 

19 

3.8 

* 

* 

10 

10 

10 

6.1 

0.9 

Lateral 

4.1 

2.6 

12.5 

10 

5.3 

2.4 

1.5 

1.4 

3.8 

5.8 

31 

3 

4 

25 

1.4 

1.2 

Vertical 

* 

4.6 

* 

1.8 

13.4 

* 

2.7 

.6 

.6 

1.1 

7.7 

0.6 

0.6 

6 

.2 

.9 

Velocity, m/sec 

















Longitudinal 

* 

3.1 

* 

1.3 

1.3 

* 

1.8 

* 

.5 

.75 

.75 

1.0 

1.0 

1.2 


.15 

Lateral 

* 

1.5 

* 

1.15 

* 

* 

.9 

* 

.5 

.65 

* 

.6 

.6 

1.2 


.30 

Vertical 

* 

.15 

* 

* 

* 

* 

.09 

.09 

.06 

* 

* 

.25 

.25 

.3 1 


.15 


Note: Asterisks indicate that no values were found. 




The actual dispersions are usually rationalized as a sum of independent random 
errors with zero means, and this sum is required to satisfy the above constraints; 
for example, the position dispersion is of the form 


6R = £ 6R. 
i 1 


and its standard deviation must not exceed the dispersion criteria, oJp, given in 
table 7.1: 


o 


6R 



(7.3) 


A traditional working model for position dispersion (e.g. , ref. 36) is that it is the 
superposition of (1) sensor errors (ground and airborne equipment), (2) navigation 
and guidance equipment errors (e.g., input data truncation, computation lags, and 
errors in computing reference trajectory and flight director or control commands), 
and (3) flight technical errors (e.g., differences between indicated and reference 
position owing to pilot or control response lags, aircraft dynamics, and external 
disturbances). A similar and commonly used rationale views the dispersions as the sum 
of independent errors arising from guidance, navigation, and control subsystems: 


6R = 6R + R + 6R n 


which are allowed to make nearly equal contributions to total dispersion. The navi- 
gation error must therefore satisfy 

°R * < V (3)1/2 0 . 4 ) 

and, similarly 

°v s »y < 3 > l/2 

Equation (7.4) is used to calculate the required navigation accuracies, which are 
included in table 7.1, except that values for VTOL landings are taken from refer- 
ence 37 and assume a 50-m square landing pad. 

The estimation accuracy of the present system is also summarized in table 7.1 for 
comparison with the required estimation accuracy. In the landing zone, the system 
meets requirements for longitudinal position (except for VTOL landings) and vertical 
position, but is deficient in lateral position. Lateral velocity requirements for 
the cases given are very nearly satisfied, but those for vertical velocity are not, 
and longitudinal velocity accuracy is midway between the VTOL and SSV requirements. 

The deficient lateral position accuracy is readily improved by relocating the MODILS 
transmitter or by replacing the MODILS navaid with the more accurate prospective 
Category III MLS system. Accuracy specifications at the runway threshold for the MLS 
are listed in table 7.1 (refs. 38 and 39) and are seen to provide sufficient lateral 
position measurement accuracy for all automatic landing operations noted in table 7.1. 
The landing zone accuracy is also compared in table 7.1 with the performance achieved 
in flight tests by the "RAINPAL" Kalman filter estimation system (ref. 3), which 
used acceleration and position sensors with significantly higher accuracy than the 
present sensors (an LTN51 inertial grade IMU and three appropriately sited Cubic range 
transponders). This system yields significantly greater accuracy than the present 
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one for the longitudinal and lateral coordinates and satisfies all automatic landing 
requirements for these axes. Finally, estimation accuracy on a 3° glide slope is 
compared with requirements for Category II automatic approaches in table 7.1. The 
known position and lateral velocity requirements are readily met. The longitudinal 
velocity requirement applies to airspeed rather than ground velocity and is met by 
the airspeed sensor, which has an accuracy of 0.6 m/sec. 

Criteria for automatically controlled area navigation operations are incompletely 
known. Position tracking accuracy criteria for area navigation are available; they 
are designed to accommodate VOR/DME and conventional altimetry and are satisfied by 
the present estimation system. However, several problems arise from the off-nominal 
and higher-frequency nominal errors of VOR/DME or TACAN which limit the usefulness of 
these navaids in support of automatic flight operations. Local experience with these 
navaids indicates that various types of error abnormalities are likely to occur on an 
approach, such as (1) extended periods of signal dropout and (2) periods of signal 
hangup or with noise standard deviation much larger than nominal. Signal dropouts 
result in dead reckoning, with linear divergence in position estimation and tracking 
errors at a rate that depends on the duration of dropout and on the IMU accuracy. 

For example, 90 sec of dropout and 500 m of drift have been experienced in the vicin- 
ity of the runway with the present sensor hardware. The possible costs of dropout 
include a reduction in safety and the maneuvering required to return to the reference 
path when navaid data are regained. Solutions for this difficulty include sensor 
changes to an inertial grade IMU or to position navaids with acceptable signal- 
dropout-duration statistics. Signal hangups pass through the present data rejection 
logic and result in potentially significant filter output transients at the end of the 
hangup, with corresponding control activity and tracking error transients. Abnormally 
noisy data are overweighted by the filter and result in higher-frequency, higher- 
amplitude velocity estimation errors and related control activity and velocity track- 
ing errors. These errors have their principal effects on ride quality and control 
activity during automatically controlled flight rather than on position tracking 
errors. In addition, the results in appendix D indicate that the nominal VOR/DME and 
TACAN accuracies result in significant longitudinal (in 4D guidance) and lateral con- 
trol activity because of estimation errors at moderate values of control bandwidth, 
but that accuracies corresponding to MODILS would reduce this activity to negligible 
levels when combined with sufficient IMU accuracy. 

The satellite-based NAVSTAR Global Positioning System (GPS), which is under 
development, offers a potential alternative to TACAN and VOR/DME (refs. 40-42). 
Projected receiver cost is comparable to VOR/DME receivers and the system accommodates 
both military and civilian users. The accuracy specification for civilian use results 
in the position fix accuracy (coarse acquisition data) given in the table below. 

Error Standard deviation, m 

Bias 15.6 

Random 15.9 

Total (RSS) 21.9 

This accuracy applies to all three axes, is independent of location and heading, and 
is equivalent to MODILS and MLS accuracy at their coverage boundaries. For the longi- 
tudinal and lateral axes this accuracy would significantly increase bandwidth capabil- 
ity and tracking accuracy; it would also reduce control activity, compared to that 
obtained here with VORTAC, and would provide excellent position data for automatic 
control globally as well as in terminal area operations. 
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8. CONCLUSIONS 


This report considers the optimization and accuracy of a Kalman filter estima- 
tion algorithm for an integrated terminal area navigation system, using sensors and 
components representative of those expected to be commonly available on aircraft with 
IFR and RNAV capability. 

The goal of the optimization was to minimize computation time requirements while 
maintaining accuracy near the maximum that could be realized from the given sensors. 

A feasible Kalman filter algorithm that satisfied this goal was obtained using various 
effective applications devices, including (1) the square root filter formulation to 
maintain computational accuracy; (2) exponentially correlated random process models 
for error states with unknown dynamics to maintain realistic measurement weighting; 

(3) partitioning the computations for measurement compression and multirate execution 
to minimize the rate of processing scalar measurements; and (4) the elimination of 
poorly observable state variables and partitioning the states and measurements for 
independent filtering to minimize the computation time for processing scalar measure- 
ments. For the measurement compression it was found that measurements could be 
accumulated at 10 Hz, compressed to a single measurement for each data type, and 
processed at 1 Hz, with only negligible loss of estimation accuracy and an order of 
magnitude reduction from the required processing rate without measurement compression. 
The time required to process scalar measurements rises with the cube of the system 
order. Elimination of poorly observable states minimizes this order, with only negli- 
gible loss of accuracy for the retained states; it resulted in a 14-state system. 
Partitioning exploits the lower execution time of independent, lower-order filters; 
it resulted in two filters associated with the horizontal-plane motion (10 states) 
and vertical motion (4 states), with a 60% reduction in execution time and only a 
negligible loss of accuracy. 

The system’s accuracy throughout the terminal area, its dependence on sensor 
accuracy, its effects on automatic trajectory tracking accuracy, and its suitability 
for terminal area operations were the subjects of extensive simulation study. 

Results of the study indicate the following conclusions. 

Estimation accuracy exhibits complex behavior in the terminal area, including 
(1) order-of-magnitude differences in position accuracy with control axis at any point 
on an approach and with location along an approach because of the range-angle nature 
of the position data; (2) significant accuracy variations in the horizontal plane 
motion states during maneuvering because of gyroscopic error transients; and 
(3) indistinguishable crosswind, sideslip, and heading measurement errors. Accuracy 
is best for the vertical axis motion, which shows simple near-stationary errors, and 
during the favorable conditions when on the glide slope and in the landing zone. 


The filter adds various useful capabilities to those of the sensors alone, 
including (1) estimation of velocity and wind, (2) improved position and acceleration 
accuracy through noise smoothing and calibration of measurement biases, and (3) veloc- 
ity error stabilization in dead reckoning. Significant further improvements do not 
appear available from further development of the filter formulation but can be 
obtained from better sensors. In this regard, position accuracy is sensitive prin- 
cipally to position sensor accuracy, velocity to both position and IMU accuracy, 
acceleration and attitude to gyroscope accuracy, and wind to velocity estimation and 
heading measurement accuracy and to lateral air velocity or sideslip sensing. 
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The operational suitability of the estimator depends on the trajectory tracking 
errors and control activity excited by the estimation errors. Analysis indicates 
that in an automatic control system these effects depend primarily on the frequency 
of the errors relative to the control bandwidth: (1) position tracking errors depend 

on translational state estimation errors at frequencies below the control bandwidth, 
(2) velocity tracking depends on estimation errors at frequencies around the control 
bandwidth, and (3) control activity depends on the estimation error spectrum at and 
above the control bandwidth. Simulation results show (1) that position tracking 
errors exceed position estimation errors significantly during large changes in posi- 
tion estimation accuracy at navaid switches and during maneuvering but are otherwise 
the same size as position estimation accuracy; (2) that velocity tracking errors 
usually exceed velocity estimation errors, except that the control converts velocity 
estimation biases to position hangoffs; and (3) that control activity tends to be 
stationary on all axes throughout the approach, is largest for the horizontal plane 
control, and is sensitive to both position navaid and IMU accuracies. 

The estimation system accuracy is satisfactory or marginal relative to safety 
criteria for various autoland operations (CTOL, STOL, SSV) , Category II CTOL 
approaches, and for CTOL RNAV based on the existing VOR/DME facilities, but is 
expected to be inadequate for future VTOL operations. Additional criteria for its 
use in fully automatic control are not available, but this report indicates the 
available estimation accuracy and its effects on the performance of an automatic (4D) 
control system for use in future evaluations. These results also indicate that com- 
prehensive and usable improvements in translational state estimation and automatic 
control for the full terminal area passenger operational envelope are potentially 
available by replacing VORTAC, accelerometers, and gyroscopes with GPS and a Schuler 
tuned IMU of sufficient accuracy, and adding air velocity vector sensing. 


Ames Research Center 

National Aeronautics and Space Administration 

Moffett Field, California 94035, October 15, 1982. 
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APPENDIX A 


POSITION FIX ACCURACY 

Aircraft position can be calculated from three simultaneous measurements of 
independent functions of position. These can be represented in this discussion as 

Y i = W + Y i ± = 1,2,3 

or, more briefly, as the measurement vector 

Y = H(R r ) + Y (Al) 

where the errors, OL) , can be assumed independent random variables with zero means. 

An estimate or position fix, R r , is obtained as the solution of equation (Al) , 
neglecting the noise; that is, satisfies 

I = H(R r ) (A2) 

The position fix error is defined by 

K ~ R r (A3) 

The relation between measurement errors and position fix errors is obtained by substi- 
tuting equation (A3) into (Al), expanding in a Taylor series in R r to first order, 
and subtracting equation (A2). This yields 

X » J(R r )R r (A4) 

where 

J = [\H] 

Here, J is the Jacobian matrix of H in which the ith row is the position gradient 
of the ith measurement type, and J is nonsingular, since it was assumed that a 
solution of equation (A2) existed. In that case, the deterministic relation between 
position-fix and measurement errors is 

R r = J" 1 ! (A5) 

and the relationship between their statistics can be obtained by applying the expec- 
tation operation, 

E[R r R Y ] = J _1 E[Y Y T ]J -T 

Or, more briefly, 

P r = J _1 QJ -T (A6) 
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where Q is diagonal, since the measurement errors C^±) are independent 


Q = diag(a i ) 

and P r is the covariance matrix of position estimation errors and contains the 
variances and correlations among x,y,z which are denoted 


P r " 


From equation (A6) , the diagonal terms of P r can be given in terms of the elements 


{a^} of J -1 * 


2 2 
J x = a ! 

j 2 = a 2 

y a 2 

j 2 = a 2 - 
Z 3 11 


£- 

X 

^xy 

q 

n XZ 

W 

a 2 

y 

q 

yz 

XZ 

^yz 


of 

Pr 

can be 


2 2 
a i2°2 

+ a i 3 

-; + 

0 2 

a 22°2 

+ a 2 3 

'5 + 

2 2 
a 32°2 

+ a 33 


,2 v 


(A7) 


Thus, the matrix J _1 (R r ) maps measurement errors Y into position errors R r . In 
general this mapping is position-dependent so that position fix accuracy varies with 
location in the terminal area. 

Analogous formulas for the path axes position coordinates can be obtained by 
recalling that 

R p - T pr <Y -*V )S r 


Equation (A4) then becomes 


X - J(E r )T pr (Y,* v )R p 


Thus, the Jacobian matrix for these coordinates is defined by 


J = J (R )T (y,ip ) 
p ' r' pr v' 


(A8) 


This matrix and the resulting position coordinate errors Rp are functions of both 
location and velocity vector direction, but since Tp r is orthogonal, the total 
position error is unchanged ( | Rp | = | R r | ) . The covariance matrix for Rp is 


P = J -1 QJ T 
P P P 


(A9) 


The matrices P r ,Pp define the position fix accuracy obtained. In the following 
analysis, equations (A4), (A6), (A8) , and (A9) are used to evaluate the accuracy 
obtained, using several combinations of measurement types of interest in the terminal 
area. 
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POSITION FIX ACCURACY WITH TACAN AND BAROALTIMETER 


This set of measurements is relied on for position information everywhere in the 
terminal area except in the region of MODILS coverage. Their position gradients were 
derived from the measurement functions listed in table 2.3 and figure 2.2 and results 
are given in table Al. The Jacobian matrix for these data types is readily assembled 
from the expressions in table Al, and its inverse is then used in equation (A5) to 
obtain the following relation between position and measurement errors: 


cos 
^T- = l Sln 


tr 

cos EL 


-sin /-cos tan EL fc 

cos I V jtb + ( - sin *t tan EL t |\b 


(A10) 


-1 


where ip t and EL t are the heading and elevation of the line of sight to the TACAN 
station. The contributions of range, bearing, and altimeter errors to position errors 
have been separated in equation (A10). As seen, the range error contribution lies 
along the line of sight projected in the ground plane, and it increases with elevation 
above the TACAN station; the bearing error effect is orthogonal to this in the ground 
plane and increases with distance from the TACAN; and the altimeter error defines z 
and contributes some error in the same direction as the range error, which increases 
with elevation above the TACAN station. 


Formulas for estimation error variances are obtained using equations (A7) 
and (A10): 


Q x = cos2 ’J't 


a 2 = sin 2 \b. 

y 


tr o 

H tan EL 

2 ttt t hb 


cos EL 


tr 

1- tan EL a/", 

2 -i-rr t hb 


cos EL 


+ (sin ip. d a., ) : 
- T t xyt tb 


+ (cos ip. d .a, )' 
r t xyt tb 


(All) 


a = a,, 
z hb 


and, further, the variance of the total ground plane position error is 


Oj = a 2 + a 2 
d x y 


^cos EL^ + ^ d xyt°tb^ + ^ tan EL t°hb^ 


The gross horizontal plane accuracy is indicated by a^, whose variation with dis- 
tance is shown in figure Al(a) for a fixed altitude (600 m) . Position accuracy 
improves from 500 m at a distance of 50 km to 75 m at a distance of 2 km. These 
errors are dominated by biases which may be greater or smaller on any given approach 
and which result in correspondingly better or poorer measurement accuracy. Results 
in figure Al for the cases that (1) all biases are small compared to the remaining 
errors and (2) all have 3 ct values indicate the large range of accuracies possible, 
depending on the bias sample values. Accuracy for V0RTAC varies from 330 m to 
1,800 m over the terminal area and is much poorer than that for TACAN. 


89 





2 



DISTANCE FROM VORTAC STATION 


(a) Ground plane position accuracy. 


ff y = 400 m 350 300 250 200 150 100 75 75 100 150 200 250 300 350 400 450 



(b) o x ,Oy contours in terminal area. 

Figure Al.- Position accuracy using TACAN and baroaltimeter measurements 
(cr tr = 70 m, = 0.58°, = 30.5 m, z = -600 m) . 
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ACCURACY ~ GROUND PLANE POSITION, m 



ACCURACY ~ RUNWAY AXIS COORDINATES, m 



ACCURACY ~ PATH AXIS COORDINATES, m 



(d) Accuracy along a STOL approach. 
Figure Al.- Concluded. 
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A terminal area map of accuracy for the x,y coordinates (fig. Al(b)) shows the 
asymptotic behavior at large distances from the station. The contour patterns for 
both coordinates are symmetric about axes through the TACAN station parallel to the 
runway x and y axes, and these two patterns are related by a 90° rotation. As seen, 
a x (<Jy) is approximately constant along lines parallel to the x-axis (y-axis) and 
reaches a minimum of a tr (70 m) on the TACAN x-axis (y-axis). Contours for a x ,a y 
in the region close to the runway (fig. Al(c)) show more complex patterns but have the 
same symmetry and rotation properties noted above. The contours for cr x ,a y - a tr 
enclose a large region around the runway within which accuracy for both x and y 
exceeds a tr . This region contains the final portion of the reference STOL approach 
shown in the figure. Position accuracy varies with location in the terminal area and 
differs for the x,y coordinates at a given location so that an approach path flown 
through these patterns can encounter large variations in the accuracy of the measured 
coordinates. For example, accuracy time histories for the STOL approach (fig. Al(d)) 
show nearly constant ground plane position accuracy, but the x,y coordinate 

accuracies show large variations associated with crossings of the TACAN y-axis (these 
excite the estimator’s transient responses), and even greater variability for the 
path axes coordinates as a result of turns (these excite the control system’s tran- 
sient response). 


POSITION ACCURACY WITH MODILS AND BAROALTIMETERS 


In this application, four measurements are used to calculate R r ; x,y are cal- 
culated from barometric altitude and MODILS range and azimuth, and then z is com- 
puted from x,y» and MODILS elevation. This solution method follows that of refer- 
ence 14 and allows computation of x,y independent of the elevation measurement, 
whose volume of coverage is interior to that of the range azimuth measurements. 

The accuracy obtained from this solution method can be computed from the position 
gradients listed in table A1 . First, the Jacobian matrix for the measurements used 
to compute x,y is 


J = 


rh n 

mr 


ma 


L h hbJ 


(A12) 


The relation of x,y to measurement errors is given from the inverse of equa- 
tion (A12) as 


GH 


-cos Az/cos EL \ / tan Az/cos EL \ / tan EL \ 

m )Y+( ) d Y + ( ) 

-sin Az / mr \ -1 / XZm ma \ 0 / 


hb 


(A13) 


where Az and EL m (aircraft elevation above the azimuth antenna) and ^ xzin are 

defined in table A1 . At the small azimuth angles for which the MODILS system is in 
use, the range error contributes principally to x, and this increases with elevation; 
azimuth error contributes principally to y, and this grows with distance from the 
antenna; and the altimeter error appears solely in x and is elevation-dependent. 

At the low elevations of normal approach paths , the altimeter error has only 
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negligible effect on the accuracy of the ground plane coordinates. The corresponding 
position error variances are obtained using equations (A7) and (A13): 


^2 _ / cos Az \ / tan Az ^ \ 

x l cos EL mrj l cos EL xzm ma / 

\ m / \ m / 


+ (tan 


E Vhb > 2 


a 2 = (sin Az a ) 2 + (d a )' 
y mr 7 xzm ma 


(A14) 


q = E[xy ] 
xy J 


(cos 2 Aza 2 - d 2 a 2 )tan Az/cos EL 
mr xzm ma m 


/ 


The error correlation q X y is included for use below in the analysis of z. 


Second, z for the solution method defined above is obtained using the MODILS 
elevation gradient in the expression 


Y 

me 


(h ,R)=ax + ay + az 
me r x y J z 


where a x , ay, and a z refer to the elevation-gradient components in table A1 . Solv- 
ing for z yields 



(A15) 


Expressions for x,y, given by equation (A13), can be substituted into equation (A15) 
to give z in terms of the independent measurement errors. A simpler computation is 
provided by forming E[z 2 ] from equation (A15) to obtain 



The distance from the elevation antenna d e is introduced to analyze the relative 
magnitudes of the terms in equation (A16); since elevation coverage is restricted to 
small elevation angles, d e a z is of the order of 1, and (a x /a z ) , (ay /a z ) << 1 in this 
expression. The accuracy therefore depends principally on elevation accuracy and 
decreases with distance from the antenna site. 

Position accuracy was computed from equations (A14) , (A16) , and (2.7), using the 
parameter values and MODILS coverage limits given in table 2.3. Accuracy for the x 
coordinate is nearly constant over the region of coverage and is given by the range 
error: 


The x-coordinate accuracy using TACAN is also nearly constant in this region and is 
given by the TACAN range error, but is much larger (70 m) so that a large jump in a x 
occurs on entering the MODILS range-azimuth coverage. Maps of Oy and a z over the 
ground plane are given in figure A2; values outside the MODILS coverage correspond to 
TACAN-baroaltimeter use. Within the MODILS region, Qy decreases an order of 
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magnitude, from 50 m at maximum distance, to 5 m in the landing zone, and this is 
everywhere better than can be achieved with TACAN. Similar behavior occurs for o z 
in the region of elevation coverage where a z decreases from 15 m to 2 m at the 
runway . 

The effects of altitude on accuracy are shown in figure A3 in which the runway 
centerplane is mapped. a z is insensitive to altitude, and Oy shows some loss of 
accuracy with altitude or glide slope angle. Outside MODILS coverage, accuracy corre- 
sponds to TACAN use; for Oy, the accuracy contours in this region are complex and 
include a vertical plane cut of the cone of confusion. 
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o z = 30.5 m 



Figure A3.- Position-fix accuracy in MODILS coverage region and vicinity 

runway center-plane map. 







APPENDIX B 


SIMULATION MODELS FOR VERTICAL AND DIRECTIONAL GYROSCOPES 


This appendix provides details of the simulation models of the vertical and 
directional gyroscopes; the models were used to generate attitude measurement errors 
in this study. The vertical gyroscope has two degrees of freedom, with its outer 
gimbal axis fixed in the aircraft parallel to the body longitudinal axis and its spin 
axis controlled to align with the local vertical. Aircraft pitch and roll angles are 
measured directly as the inner and outer gimbal angles, respectively, and contain 
errors determined by the misalignment of the spin axis from the local vertical. The 
misalignment is governed by a dynamic equation and depends on aircraft acceleration 
and attitude histories. The directional gyroscope has two degrees of freedom, with 
its outer gimbal axis parallel to the body-normal axis; its spin axis is controlled 
to align with magnetic north. Heading is measured as the outer gimbal angle and 
contains errors determined by the aircraft attitude and misalignment of the spin axis 
from magnetic north. These errors are also deterministic functions of the aircraft 
motion. 


VERTICAL GYROSCOPE 


Axis Frames, Orientation Angles, and Transformations 

Earth, local vertical, and body-axes systems are the usual orthogonal reference 

frames; they are illustrated in figures B1 and B2. These frames are denoted F e , F v , 

and Ffo and are defined by the unit vectors 

F e : 

V {n,e,d} 

F b : 

The transformation from local vertical to body axes coordinates is given by 

cos 0 cos ip cos 0 sin ip -sin 0 " 

sin <J> sin 0 cos ^ - cos <j> sin ip sin <J> sin 0 sin ^+cos < p cos ip sin < p cos 0 

cos (}) sin 0 cos ijj+sin <p sin ^ cos (p sin 0 sin ip - sin <j> cos cos <p cos 0 

(Bl) 

The vertical gyroscope is mounted in the aircraft with its outer gimbal axis 
along the body longitudinal axis, i^. A gyroscope axis frame is illustrated in 
figure B3 and is defined as 



F 

S 


{£ x >iL i >iL 3 } 
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F e = Earth axes; Earth centered and rotating with Earth, referenced to 

equatorial plane and Greenwich meridian 

F v = {N,JE,D^} local vertical axes; north, east, down, with origin at the gyroscope 
location given by spherical coordinates (R,£,X) 

R Earth radius 

£ longitude angle 

A latitude angle 

a) Earth’s angular velocity 

Figure Bl.- Earth and local vertical reference axes. 


where ^ and _gg are the inner gimbal axis and the rotor spin axis of the gyroscope, 
and fLx i- s orthogonal to the two vectors. This frame is obtained from F b by suc- 
cessive rotations through the outer gimbal angle, G Q , counterclockwise about jL b , and 
through the inner gimbal angle, G^, counterclockwise about _g^. The transformation 
from F b to F g is, therefore. 


T gb " M-<VM- G o> ‘ 


cos G. sin G. sin G 
1 10 


cos G 


sin G. cos G. sin G 
1 10 


sin G. cos G 
1 o 


-sin G 


cos G. cos G 
1 o 


The spin axis is located in the local vertical frame by angles E and F, shown 
in figure B3. These angles measure the misalignment of the spin axis from the local 
vertical. It is useful to construct the misalignment axis frame, F m , as 


F 

m 


= <Vin’-£s } 


which is oriented relative to F v by angles E and F; from the geometry shown in 
figure B3, the transformation from F v to F m is given by 
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= {i fe , jfa } axis system fixed to the aircraft 

V* heading; angle from north in the horizontal plane 

6 pitch; angle between an d the horizontal plane 

^ roll; angle of rotation about from the horizontal plane to 

Figure B2.- Aircraft body axes and Euler angles. 


cos F sin F sin E -sin F cos E" 

T mv = E 2 (F)E 1 (E) = 0 cos E sin E (B3) 

_sin F -sin E cos F cos E cos F 

Finally, the gyroscope frame, Fg, can be obtained by a rotation about the spin axis. 
The required angle is denoted G s ; then, 

"cos G s sin G s Cf 

T gm = E 3 (G s ) = “ sin G s cos G s 0 (B4) 

0 0 1 


It is noted here for later use in vector operations that the elements of a 
transformation between orthogonal reference frames, say from {a 1 ,a 2 ,a 3 } to 
{b a ,b 2 ,b 3 }, are the inner products of the unit vectors 



(B5) 


and that the ith row of the T a b expresses the vector aj_ as a linear combination 
of the axis vectors of 


101 



F g = jSx'Si'SsJ GYRO AXES, g j( g $ ARE THE 

INNER GIMBAL AND GYRO SPIN 
AXES, RESPECTIVELY 

G OUTER GIMBAL ROTATION 
ANGLE (ABOUT i b » 

G- INNER GIMBAL ROTATION 
‘ ANGLE (ABOUT gj) 


(a) Orientation in aircraft body axes. 

F m “ |!m'jm'5s| 

E MISALIGNMENT AXES 
F ROTATION ANGLE ABOUT N 
Gs ROTATION ANGLE ABOUT j m 
ROTATION ANGLE ABOUT g s 



(b) Orientation in local vertical reference frame. 

Figure B3.- Vertical gyroscope: axes and orientation angles. 
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(B6) 


— i = + ^i’i.2^2 + (—± ’—3 ^3 

Similarly, the jth column of T a b gives bj as a linear combination of 

{a^a^aj}. 


Gimbal Angles 

The angles G 0 , G^, and G s can be expressed in terms of the aircraft attitude 
angles and spin axis misalignment angles as follows. The transformation T gv can be 
formed from alternative sets of angles, 

T gv ■ T gm (G S )T mv <E - F ) - T gb (G i- G o )T bv«' e '« < B7 > 

and this can be rearranged as 

T gb (G i- G o>V G s> ■ T bvM’ 6 '« T Iv (E ’ F > 

or 

E 1 (G o )E 2 (G i )E 3 (G s ) = E 1 (<fr)E 2 (0)E 3 (#)E 1 (-E)E 2 (-F) (B8) 

From equation (B8) , {G 0 ,G-^,G S } are identical to {cf>,0,^}, respectively, if the spin 
axis is aligned with the local vertical (E = F =0), and conversely. Solutions are 
readily obtained after expanding the left-hand side of equation (B8) and using ele- 
ments of the right-hand side as needed; in an abbreviated form these are 

tan G s - <i b ,j 1n >/<l b ,i n > 

sin G ± = (B9) 

tan G 0 = <J b .£ 8 >/<k b .£ 8 > 

The inner products are readily expanded in terms of <j), 6, 4s E, and F, using the rows 
°f ^bv> Tmv ( ec is. (Bl) and (B3)) and equations (B5) and (B6) . The results are listed 
in the simulation equation summary (see table Bl). 

Misalignment Equation 

Equations governing the spin axis misalignment angles can be obtained from the 
equations governing the forced gyroscope dynamics. These satisfy 

R = “g/! ® £s (BIO) 

where 

= angular velocity of Fg relative to inertial space 
£ = torques applied to the gyroscope 
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Figure B4.- Simulation of vertical gyroscope: computational flow diagram 





TABLE Bl.- SIMULATION OF VERTICAL GYROSCOPE: EQUATION SUMMARY 


Total acceleration and spin-axis misalignment from apparent vertical 


fg cos G cos F sin G cos E + cos G o sin F sin E sin G sin E - cos G sin F cos £ 

* 5> s s s s 

fg = fg 2 = -sin G s cos F cos G q cos E - sin G q sin F sin E cos G g sin E + sin G g sin F cos E 


-cos F sin E 


cos F cos E 


a 0 “ f 8i/| £ g 3 ! 
“i “ f g 2 / I f S 3 I 


Gimbal angle control rates, rad/sec 
P i = P E (0t i ) + 1 * 92xl0_5 


P o = P E (cx 0 ) + 3.84x10' 


-5.7596x10' a 


a < 0.017453 


P E = <5.7596x10 sign (a) ja| e [0.017453, c], (c - 0.1047 for V ± , c = 0.0436 for P Q ) 


Spin-axis precession relative to local vertical 


fl-E 


-cos G /cos F -sin G cos G 
s s 


cos G„ cos 
s 


cos F -(cos >. 

+ to 

G. p n 

i J L°J L 


+ sin X cos E tan F)" 


sin X sin E 


^l-sin ip(l + tan X cos E tan F) + cos sin E tan F 
P I sin i|> tan X sin E + cos tji cos E 

Gimbal angles and gyroscope pitch and roll outputs 

sin G^ = cos F sin E sin ^ cos 6 + cos F cos E sin 0 - sin F cos ij; cos 6 E ± y 


_ cosE cos 8 sin $ - sin E(sin ^ sin 9 sin $ + cos if; cos tft) + tan F(cos \\> sin 6s in - sin^ cos <f> ) 


° cos E cos 6 cos <{> - sin E(sin ij; sin 0 cos <J) - cos ^ sin <J>) + tanF(cos ^ sin 0c os <j>+ sini|) sin <}>) 


t = G 

o 


9 g = G i 


0 = 0 - 0 
g 


cos E sin cos 0 + sin E sin 0 


s cos F cos cos 0 + sin F(sin E sin ^ cos 0 + cos E sin 0) 


Initialization 


sin F - -(sin G^ cos 0 - cos cos G^ sin 0)cos ^ - sin cos G^ sin iJj 
( sin G. cos 0 - cos 5 cos G. sin 0)sin ip - sin $ cos G. cos ^ 

ctn V - ± ° 1 O 1 


FS! 2 


Eei I 


cos <f> sin if; + sin $ sin 0 cos 


(cos G. cos 9 + cos sin G. sin 0)cos if> - sin $ sin G, sin 


© 



The torque is expressed as an angular rate (torque per unit angular momentum) in 
equation (BIO) and throughout this appendix. From the geometry of figures B1-B3, the 
angular velocity can be given as the sum of independent rotational rates, 

= (a) + £)K - X E + EN + (Bll) 

where £,X are given from the aircraft velocity as 

£ = V sin ip/R cos X 
X = V cos ip/R 

All unit vectors in equation (BID can be given in terms of the unit vectors defining 
F m , using T mv (eqs . (B3) and (B6)) as needed, 

N = cos F + sin F 

E = sin F sin E + cos E - sin E cos F 
K = cos X N - sin X I) 

= (cos X cos F + sin X sin F cos E)i^ - sin X sin F 
+ (cos X sin F - sin X cos E cos F)_g 


and then the cross-products required in equation (BIO) are 


N ® = -cos F 


_E ® = cos E i^ - sin F sin E 


K ® _g s = -sin X sin E i^ - (cos X cos F + sin X sin F cos E)j^ 


Finally, substituting equations (BID and (B12) into (BIO) and taking inner products 
with the orthogonal vectors ijn» jjn yields the following equations for the misalign- 
ment angles: 


<E'V 

cos F 


- oo[cos X + sin X tan F cos E] 


+ £ [-sin xp ( 1 + tan X tan F cos E) + cos ip tan F sin E] 


<p,i ) + (o sin X sin E + (sin ip tan X sin E + cos \p cos E) 
— m K 


These equations are integrated in the simulation to obtain E and F, once the control 
torque £ is given. 
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Gyroscope Control Torques 

Torque is applied along the inner and outer gimbal axes, and can be written as 


£ - P& + Poh 

Noting that can be represented as 

i, = cos G.g - sin G,g 
— b i— x l-^s 

then the inner products required in equation (Bll) can be written as 


{p.i^) = -p 1 sin G g + p Q cos G i cos G g 


= Pi cos G s + P G cos G i sin g £ 


(B14) 


(B15) 


For both gimbal axes, the applied torque is due to control torques and fixed 
drifts associated with mass unbalance. These are denoted as 


and suitable drift values are 


P i ■ P iC + P iD 


P o P oC + P oD ) 


P iD = 0.0011°/sec 
P oD = 0.00225°/sec 


(B16) 


(B17) 


The control torques are functions of the misalignment of the spin axis of the 
gyroscope from the apparent vertical; the misalignment angles are sensed by orthog- 
onally located bubble levels on the rotor case and are given by 


% = <i*8 x >/( i’S s > ' 


(B18) 


where _f is the specific force acting on the gyroscope, 

f. = a - £ 

and its direction is the apparent vertical. The gyroscope axes coordinates of f_ 
needed to evaluate equation (B18) can be computed from its local vertical axes compo- 
nents, which are available from the aircraft motion simulation, using equations (B3) 
and (B4) : 


f g = E 3 (G s )E 2 ( F )E i (E)f v (B19) 

The misalignment a Q is controlled by precessing the inner gimbal using torque along 
the outer gimbal axis in accordance with equation (BIO). Similarly, <Xj[ is controlled 


107 




by precessing the outer gimbal, using torque 
along The gimbal control law for both 

gimbal axes is modeled approximately by the 
functional form shown in sketch I. 




' -(b/a)a 
- -b sign(a) 
0 



< a 

€ [a,c] 


o 


> c 


(B20) 


where the parameter values are 

a i = = i° 

b-£ = b Q = 0.033° /sec = maximum gimbal rates 
c G = 2.5° 


c 


i 


= 6 ° 


Sketch I The control is linear, with saturation at a 

maximum gimbal rate of 0.033° /sec and with 
disengagement at large misalignments from the apparent vertical. In static equilib- 
rium flight, the apparent and local vertical are aligned and the misalignments a Q and 
remain small and within their cutoff values. During aircraft maneuvering, sub- 
stantial angles between the two verticals develop, but they develop at a much higher 
rate (by one to two orders of magnitude) than the maximum gimbal control rates. 

Thus, for accelerations sufficiently large to exceed the cutoff misalignment (0.05 g 
longitudinally and 0.1 g laterally) the gyroscope control is disengaged before sig- 
nificant misalignments from the local vertical, E and F, can be developed by the con- 
trol torques. After disengagement, the gyroscope is uncontrolled about the disengaged 
axis, but the rate of precession is sufficiently small that very little error develops 
during typical maneuver durations. 


Simulation Summary 

The vertical gyroscope can be simulated using equations (B9), (B13), (B16), 
and (B20), with aircraft accelerations, Euler attitude angles, and latitude as 
inputs. A computational flow chart is presented in figure B4 and an equation summary 
for the simulation is given in table Bl. The gyroscope can be initialized at nonzero 
errors, <j>,0, by calculating appropriate initial values for the gimbal and misalign- 
ment angles. Formulas for these values were derived from equation (B8) and are 
included in the equation summary. 


DIRECTIONAL GYROSCOPE 


Axes and Transformations 

The directional gyroscope is mounted in the aircraft with the outer gimbal axis 
fixed along the body normal axis, k^. A gyroscopic axis frame is illustrated in 
figure B5 and is defined as 


F g = >£i’£ z } 
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F g " {is' Si' ?z} GYR0 AXES F0R 

dirIctional gyro 

G- OUTER GIMBAL ROTATION ANGLE 
ABOUT k b 

G: INNER GIMBAL ROTATION ANGLE 
ABOUT gj 


F m = {is* im' hm) DIRECTIONAL GYRO 
MISALIGNMENT AXES 

A SPIN AXIS AZIMUTH FROM NORTH 
ABOUT D 

T SPIN AXIS TILT ANGLE FROM 
HORIZONTAL PLANE 

G ROTATION FROM HORIZONTAL PLANE 
TO gj ABOUT g s 




(a) Orientation in body axes* 


(b) Orientation in local vertical 
reference frame. 


Figure B5.- Directional gyroscope: axes and orientation angles. 


where _gg and are the rotor spin and inner gimbal axes and is orthogonal to 

these two vectors. This frame is obtained from the body— axes frame, F^, by successive 
counterclockwise rotations through G 0 about and G ± about ^ (see fig. B5); 

hence 


T 


gb 


e 2 (-g.)e 3 (-g o ) 


sin cos 
sin G 

o 

-sin G. cos 
1 


G 

o 


G 

o 


-cos G. sin G 
1 o 

cos G 

o 

sin G. sin G 
1 o 


sin G i 
0 


cos G. 

X 


(B21) 


The spin axis is located in the local vertical frame by the azimuth (A) and 
tilt (T) angles shown in figure B5. These angles measure the misalignment of the spin 
axis from local magnetic north. A misalignment axis frame is defined as 
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F m ' <VW 


and is oriented relative to F v by the successive rotations A about D and T 
about hence 


T 

mv 


= e 2 (t)e 3 (a) = 


cos T cos A 
-sin A 
sin T cos A 


cos T sin A 
cos A 

sin T sin A 


-sin T 

0 

cos T 


(B22) 


Finally, the gyroscopic frame, F g , can be obtained from F m by a rotation G s about 
the spin axis which locates ^ in the misalignment frame: 


no 0 


T 

gm 


MG S > 


0 

COS 

G 

s 

sin 

G 

s 

0 

-sin 

G 

cos 

G 


(B23) 


Gimbal Angles 


Expressions for G Q , Gj_, and G s in terms of aircraft attitudes and spin axis 
misalignment angles are obtained from the transformation identity 


T 

gv 


= T (G ) T (A,T) 
gm s mv 


- VV G o )T bv ( 'M’' ,,) 


which can be rearranged as 

E, (G )E 2 (G.)E 1 (G ) = E 1 (^)E 2 (e)E 3 (4 ) )E 3 (-A)E 2 (-T) 

O X s 


(B24) 


Solutions for G 0 , G-^, and G s are given in the equation summary, table B2. In 
abbreviated form, these are 


tan G o = -<i b .^ s >/<i b ,l s > 

sin G ± = (B25) 

tan G s = 

The inner products can be expanded using equations (Bl) and (B20). For G 0 this 
gives 


tan G 

o 


cos (j) sin(^ - A) - sin cj) sin 6 cos(^ - A) - tan T sin <fr cos 6 
cos 0 cos(ty - A) + tan T sin 0 


(B26) 


In the case in which <J> and 0 are zero, the measurement error depends only on the 
azimuth misalignment G 0 = i|> - A and is zero if A is nulled. For flight with 
* = 0, 


110 



TABLE B2 . - SIMULATION OF DIRECTIONAL GYROSCOPE: EQUATION SUMMARY 


Heading and heading misalignment angles 


cos T cos A 


cos T sin A 


fg “ ^ j* = sin sin T cos -A - cos G„ sin A sin G sin T sin A + cos G cos A sin G cos T f 

S S S S S 

f e cos G sin T cos A + sin G„ sin A cos G„ sin T sin A - sin G cos A cos G cos T 


- U! Uo 


cos 0[sin 0 (cos 0 - 1) - tan sin 0 
tan ^c 1 + sin 0[sin 0(cos 0 - 1) - tan y^ sin 


c . - + if/ - G 

0 Y Y c o 

Gimbal angle control rates, rad/sec 

p i = p c (e * ) + 3 - wlx10 ' 5 


P o = p c( a o> + 4.887X10- 5 + P 0Df 


P c (a) = 


(a,b) = 


(b/a)a [ a | i a 

b sign(a) |a| > a 

(8.73xl0 -2 , -1.92xlCT 3 rad/sec) for P 

( 1 . 745* 10“ 2 » 8.73xl0~ 4 rad/sec) for P_ 


0 | G 0 | < 1 . 745* 10 -5 rad/sec 

3.49 x 10 -4 sign(G 0 ) |G 0 | > 1.745*1CT 5 rad/sec 


Spin axis precession relative to magnetic north 


A cos 
J -£ 


G g /cos T -cos G^ sin G g /cos T I I Isin X - cos X tan T cos A 


-cos G. cos G 
i s 


i] R 1 

oJ l 


cos X sin A 


y sin 0(tan X - tan T cos A) + cos 0 tan T sin A 
sin 0 sin A + cos 0 cos A 

Gimbal angles and heading measurement 

cos sin(0 - A) - sin 0 sin 8 cos(0 - A) - tan T sin j) cos 6 
an o cos 0 cos(0 - A) + tan T sin 0 

sin G^ = cos T cos A(cos 0 sin 0 cos 0 + sin 0 sin 0) + cos T sin A(cos 0 sin 0 sin 0 - sin 0 cos 0) - sin T cos 0 cos 0 

n £ _ -[-sin A(cos 0 sin 0 cos 0 + sin 0 sin 0) + cos A(cos 0 sin 0 sin 0 - sin 0 cos 0)] 

s sin T cos A(cos 0 sin 0 cos 0 + sin 0 sin 0) + sin T sin A(cos 0 sin 0 sin 0 - sin 0 cos 0) + cos T cos 0 cos 0 


*g " *g - ♦ 


Initialization with error 0 Q and gimbal angle 

G n - 0 + 0 


tan (A - 0) = 


-cos G,^ sin G q cos 0 - sin G^ sin 0 
cos G^(cos G q cos 0 - sin G q sin 0 sin 0) + sin G^ cos 0 sin 0 


sin T = cos G^cos G q sin 0 + sin G q cos 0 sin 0) - sin G^ cos 0 cos 0 T 6 i ^ 
-sin G^ sin 0 + cos G cos 0 sin 0 

£ an Q _ 2 2 

an s sin G^(cos G q sin 0 + sin G q cos 0 sin 0) + cos G^ cos 0 cos 0 


in 


<f>, 'll X4>,v <P, o, ^ 



Figure B6.- Simulation of directional gyroscope: computational flow diagram. 




tan G 

o 


’ sinpfr - A) I 1 

cos (\p - A) + tan T tan ej cos 0 


Perfect azimuth alignment does not null the measurement error in this case, but the 
error is small and insensitive to small values of both 0 and T. For nonzero <}>, 
errors are dynamic and dominated by inaccurate control as discussed later. 


Misalignment Equations 

Equations governing the spin axis misalignment angles are again obtained from 
the dynamics of the forced gyroscope. These satisfy 


£ = W g/X ® (BIO) 

where, for the directional gyroscope, 

-g/I = (“ + - A E + AD + (B27) 

Express all unit vectors in equation (B27) in terms of F m using equation (B22) , 

N. = cos T cos A — sin A + sin T cos A 

_E = cos T sin A + cos A + sin T sin A 

D = -sin T g + cos T k 
— -“s — m 

K = cos X N — sin X I) = (cos X cos T cos A + sin X sin T) g_ 


- cos X sin A + cos X(sin T sin A - sin X cos T)!^ 


and then the cross-products required in equation (BIO) 
N ® .gg = sin T cos A + sin A k^ 

® £s = sin T sin A + cos A k^ 


are 


(B28) 


D ® jg = cos T ^ 

K ® ^ = (cos X sin T cos A - sin X cos + cos X sin A k^ 

Substitution of equations (B27) and (B28) into (BIO) , and formation of inner products 
with jm an< i iSm> yields the following equations for the misalignment angles: 
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^£.»in > 

A = - — b co(sin X - cos X tan T cos A) 

cos T 


+ [sin (tan X - tan T cos A) + cos tan T sin A] 
R 


1 


T = + co cos X sin A - [sin i|i sin A + cos cos A] ) 


(B29) 


Gyroscope Control Torques 

Torques can be applied along the inner and outer gimbal axes. 


£ - PA + Poib 


Noting that 


h = _sin G i^s + cos G ±&z 

then the inner products required in equation (B29) can be written, using equa- 
tion (B23) , as 


(P,i ) = p. cos G - p cos G, sin G 
^ s r o 1 s 


(B30) 


(p,k > = p. sin G + p cos G. cos G^ I 
— m s *o l s / 

The applied torques along both axes are due to control torques and drifts and are 
denoted as 


Pi = PiC + P iD 


P o P oC + P oD ) 


(B 16 ) 


where suitable drift values are 

p^ D = 0.002° /sec 

(° 

p = 0.0028 + { 

° U 10.02 sign(G 0 ) 


| G 0 1 < 0.001°/sec 
| G 0 | > 0.001°/sec 


(B31) 


The outer gimbal drift contains a fixed drift owing to mass unbalance and a friction 
drift owing to gimbal rate which is fixed in magnitude for rates above a threshold 
rate. The rate G Q can be approximated as the aircraft yaw angular velocity in the 
simulation. 

The outer gimbal control (see sketch J) acts to maintain the spin axis in the 
apparent horizontal plane by precessing the inner gimbal to null the misalignment 
angle , 
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a o 5 <i»^ s >/|<I»g z > l (B32) 

using the control 


b/a a Q 

Q 

o 

A 

b sign(a o ) 

a > a 


Suitable values for the saturation angle 
and maximum gimbal precession are 



Sketch J 


a = 1° 

b = 0.05°/sec 

The inner gimbal control torque acts to null the error between the outer gimbal 
angle and the apparent heading, i pQ, sensed by a pendulous magnetic sensor and is 
modeled by 


where 


f -0.022 ^ 

[-0.11 sign(e^) 


e. < 5‘ 


e. > 5‘ 

ipl 


(B33) 



The principal error in the apparent heading 

*c = * + k 


occurs in turning flight and results from detection of the direction of the magnetic 
field in the apparent horizontal plane rather than the true horizontal plane. To 
model this error, let h be the direction of Earth's magnetic field. Neglecting 
field anomalies, then 


h = cos Y D N + sin D 

where Yd is the dip angle and varies with geographic location. At Ames Research 
Center, the dip angle is 61.8°. Assume that both the sensor pendulum axis and body- 
normal axis kjj align with the apparent vertical and that 0 is zero. Then the 
sensed direction of north is located in the plane of (h^k^) and is perpendicular to 
kjj, that is, along 

n = h - (h-k^) k^ 

The apparent heading angle is given from 

cos \p c = <i b ,n)/|n| = <i b >h)/|n| = cos y d cos ij//|n| 

s ^ n = - ^ib’— ^ I— I = 3b 1— I = ip cos <p cos y^ - sin <f> sin Yp)/| n | 


115 



The error angle is obtained by forming sin(^ c - yp) 9 cos(^ c - ij/) and then 

sin(^ - cos ^[sin ij;(cos <P - 1) “ tan y D sin cf> ] 

tan cos(\p^ - yp) 1 + sin ^[sin yp(cos <f>-l)-tan y d sin<J>] (^34) 

This error is illustrated in figure B7 ; it is zero at zero roll angle and otherwise 
varies sinusoidally with heading during a turn and has large amplitude, even at modest 
roll angles (10°). During turns, this error dominates the heading error used 

for the inner-gimbal control (eq. (B33)), but the control saturation limit 
(0.11°/sec), along with the sinusoidal nature of e^, limits the resulting measure- 
ment error amplitudes during typical turns at about 5° (see fig, 2.2(c) of the text). 



HEADING ~ deg 

Figure B7.- Directional gyroscope: error in sensing north during turns 

(y D = 61.8°). 


Simulation Summary 

The directional gyroscope can be simulated using equations (B25) and (B29) 
through (B34) , with aircraft accelerations and speeds, Euler attitudes, and latitude 
as inputs. A computational flow chart and equation summary for the simulation is 
given in figure B6 and table B2. 
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APPENDIX C 


ERRORS IN MEASURING RUNWAY-REFERENCED ACCELERATIONS 


The measured runway axis components of acceleration are derived from specific 
force and attitude measurements provided by body mounted accelerometers and attitude 
gyroscopes: 


am 


“r - T rb ( V VV fm b + s r 
This differs from the actual acceleration by the measurement error 


am = a + a 
r r r 


(Cl) 


(C2) 


where a r can be expressed in terms of the sensor errors by the variational equation. 


[ 3am Sam 

r r 

. , T~Z % 

3<j> 36 


3 am 


dip 


and 


+ T,f 


rb tm b 


(C3) 


3 am 

ST ^ 

r 

rb 

3(f> 

3<f) 

Sam 

ST , 

r 

rb 

30 

30 

Sam 

3T , 

r 

rb 


3T 


fm b = 


T t fm 


T, fm 


rb 


~ , T, f m 
3<p br r 


dip dip br r 

Recalling the transformation T^ r (table 2.2), 

"l 0 0 ' 

0 cos <p sin <j> 

_0 -sin <f> cos <J> 


(C4) 


T br = E 1 (cf ) )E 2 (0)E 3 (cp) = 


cos 6 0 -sin 0 

0 1 0 
sin 6 0 cos 6 


cos ip sin ip 0 
-sin ip cos ip 0 
0 0 1 


(C5) 


then the required derivatives in equation (C4) are found to be 
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3am T _ 3E 1 (<()) 

-jf- = M*> E 2<e> — af - E 1 (<D)E 2 (e)E 3 Wfm r 


sin 0 fm + sin \p cos 0 fm 

r 2 r 3 


-sin 0 fm - cos tfi cos 0 fin 
r , 


,cos 0(-sin \p fm + cos ip f m )y 

\ r i r 2^ 


cos ip fm 


8 am 


3E T (0) 


- 3 Q 1 = Eg 010 -gQ— E 2 (0)E 3 (^)fm r = sin ip fm^ 


-cos ip fm - sin ip fm 

* r i r 2 > 


f- fm 


Bam 3E, 0(0 


3ip 


8iJj 


E 3 (^)fm r = fm r 


and the complete expression for the measurement error in equation (C3) is 


am = 
r 


sin 0 fm + sin ip cos 0 fm 
r T r 

r 2 L 3 


-sin 0 fm - cos ip cos 0 fm 
r r r 

l r 3 


cos \p fm 


sin \p fm 


-fm 


fm 


L cos 0(-sin \p fm + cos \p fm ) -cos ip fm - sin ip fm 

r i r 2 r i r 2 




g 


+ T ,f 


rb tm b 


(C6) 


The first term in this expression provides the model of the dependence of acceleration 
errors on the errors of the attitude gyroscope that is used in the 15-state filter. 

It is also of interest to express the error term of the attitude gyroscopes in 
equation (C6) in level heading axes and to introduce aircraft accelerations. One way 
to do this is to express f as 


f = a - £ = ai i L + a 2 i L + (a 3 - g)!^ 


(C7) 


where {a^} are the level heading coordinates of acceleration, and to note the vector 
relations between the level heading axes and runway axes: 


i L = cos \p i r + sin \p j_ r ^ 
1 L = -sin \p i_ r + cos \p j_ r 


(C8) 
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Using (C7) and (C8) , we have 



<f,i r > 




a 1 cos ip - a 2 sin ip 
a x sin ip + a 2 cos ip 
a 3 - g 


(C9) 


The columns of the attitude error coefficient matrix in equation (C6) can now be 
given in terms of { 11 >j_L »iir ^ usin 8 equations (C7) to (C9). For example, the third 
column gives the effect of ipg, and it can be rewritten as 


am 


($ e ) 


= (-f 



+ 



g 


= [a-^-sin ip ± r + cos ip j_ r ) + a 2 (-cos ip ± r - sin ip j. r )]^ g 

= (a il L " a 2i. L ^g 


Proceed similarly for the remaining columns, except to assume 0 is a small angle, 

I a i I » I a 3 I <<: and retain only the leading terms. The resulting expression for the 

gyroscopic error term is then 


am (syro) = (gi L + a 2 k.H - (gi. L + a^S + (a^ - a 2 i L 




(CIO) 
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APPENDIX D 


EFFECT OF ESTIMATION ERRORS ON TRANSLATIONAL CONTROL 


The effect of estimation errors on control of the aircraft translational motion 
is examined in this appendix* For this purpose a simplified generic model of the 
closed loop, automatic control system is derived and used to relate trajectory track- 
ing errors and control activity to estimation errors and to generate both analytical 
results (transfer functions and frequency responses) and simulation results for these 
relationships. The simplified model represents solely the effects of translational 
state estimation errors on translational control and does not require representation 
of any details of the aircraft force and moment generation processes. The effects of 
estimation errors in other states required by the complete system (e.g., those 
required to control attitude, engine, configuration, and actuator dynamics) are 
independent . 


A SIMPLIFIED MODEL OF AUTOMATIC TRANSLATIONAL CONTROL 


A model of the complete reference trajectory tracking system is shown in fig- 
ure Dl(a). Only the translational control will be considered and, in general, this 
can be formulated as an inverse of the aircraft translational motion. Referring to 
figure Dl(a), the aircraft translational motion model begins with the control u 
which generates the accelerations F(z,u), whose integrals give the states R r ,V r . 

The control u includes those independent variables of the acceleration generation 
process which the control laws have at their disposal (e.g., for CTOL aircraft these 
are flap, engine rpm, roll angle, angle of attack, and sideslip angle) and z refers 
to the remaining variables affecting acceleration (e.g., dynamic pressure, atmospheric 
parameters, and aircraft weight). The control inverts this logical flow. It begins 
with a reference trajectory command, [ RC r (t ) , VC r (t ) , ac r ( t ) , t Q < t < tf], which is 
assumed to satisfy Newtonian mechanics, aircraft acceleration generation limits, and 
operational constraints. In the absence of disturbances, this can be executed using a 
partial inverse of the acceleration generation process, G(z,a), to calculate those 
control settings, uc(t), that generate ac r (t), at the estimated current flight 
conditions, z(t). 

Trajectory tracking errors develop from various external and internal distur- 
bances (wind turbulence, estimation errors, control errors) and their estimates, 
6Rr,6V r are fed back to obtain stable tracking. The regulator output is a corrective 
acceleration command, Aa r , which is executed by summing with ac r and passing the 
total, act r , through the plant inverse, and which imposes transient response dynamics 
on the translational tracking errors, as specified in the feedback control law. An 
appropriate control law is given in figure Dl(b); the selected dynamics (1) are 
decoupled in path axes; (2) are second order in the linear domain, with damping and 
natural frequency determined by the gains in accordance with the relations given in 
figure Dl(b); and (3) have saturation limits on the control authority Aa^^ and on 
the velocity excursions that can be used to null large position errors, Path 

axes decoupling is used in view of the nearly decoupled control over accelerations 
along these axes provided by the conventional controls governing engine rpm, roll 
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Figure Dl.- Automatic reference trajectory-tracking system: 


translational degrees of freedom. 











^va' Ta 


LINEAR DOMAIN CONTROL 

PARAMETER VALUES 


A A 


ii 

a 

CO 

<1 

r SR p tK v 6V p 

fj = 0.707, i = 1, 2, 3 

K r = 

VJ 

1^;} = {0.125,0.125,0.3} rps 

K v = 


{ 5V li mi } = {5,5,2.5}mps 



=0.1 9, i = 1,2,3 


(b) Trajectory regulator. 



(c) Acceleration response error compensation. 
Figure Dl.- Continued. 
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f 

ACCELERATION 
CONTROL ERRORS 


(d) Simplified estimation and control system model. 
Figure Dl.- Concluded. 


angle, and angle of attack. 3 Suitable control law parameter values for passenger 
operations are listed in figure D1 (b) • The natural frequency u) n determines system 
bandwidth and is a useful gain-adjusting parameter to modify system dispersions 
caused by disturbances or to maintain acceptable control activity. Its selected 
values for the present discussion are representative of passenger operations where 
w n ranges approximately from 0.05 to 0.5 rps, depending on trajectory segment, con- 
trol axis, and disturbance levels. 

The system described so far will reproduce the reference trajectory and achieve 
the specified transient response in the absence of estimation and trajectory synthesis 
errors, provided the total acceleration command act is reproduced exactly by the 
combined plant inverse, control response, and aircraft acceleration generation 


3 

However, some coupling and alteration of the path axis transient dynamics is 
introduced by Coriolis accelerations in the kinematic relations: 


6R 

= 6V - 

(w 

® 6R) + HOT's 

P 

P 


P 

6V 

= Aa 

(w 

® 6V) + HOT's 

P 

P 


P 


Nonnegligible effects occur principally during turns in combination with large hori- 
zontal plane tracking errors. 
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process; that is, provided this combination is unity. In general, this combination 
differs from unity because of identification and approximation errors in the plant 
inverse, G, and the control response dynamics. However, the plant inverse errors are 
assumed to be slowly varying compared with the bandwidth of the acceleration commands 
and can be estimated and compensated automatically by feedback of the integrated 
acceleration error; this feedback and its control law are shown in figures Dl(a) 
and Dl(c). In addition, the control response dynamics in figure Dl(a) refer to the 
remaining aircraft degrees of freedom (attitude, engine, and actuator dynamics) 
together with their appropriate control laws, sensors, and estimation logic; in pas- 
senger operations, it can be assumed that the dynamics for these subsystems are an 
order of magnitude faster than the bandwidth of the acceleration commands so that they 
can be neglected in modeling the translational motion and its control. 

The model of the translational control system can now be simplified to that shown 
in figure Dl(d). Errors in controlling acceleration that result from control 
response dynamics, plant inversion errors, errors in estimating the variables of the 
inner loop control, and external disturbances are collectively represented by the 
additive disturbance, f. However, f can be ignored here since its effects on 
trajectory control are independent of those owing to translational state estimation 
errors, and then the simulation is easily expanded to include the control system of 
figures Dl(b), Dl(c), and Dl(d) in order to study these latter effects. 

The design structure of the automatic control, which is given in figure Dl(a) 
and which leads to the simplified model in figure Dl(d), has been proposed and flight 
tested at Ames in recent years (refs. 43 and 44). Conventional automatic control 
systems based on a linear perturbation model associated with a nominal operating 
point can be rearranged to give the structure shown in figure Dl(d) for the neighbor- 
hood of its operating point so that the results given next for relations between 
estimation errors and system dispersions apply to conventional designs as well. 


TRACKING ERRORS AND CONTROL ACTIVITY CAUSED BY ESTIMATION ERRORS 


Transfer functions relating system output variables to estimation errors are 
readily derived from the model in figures Dl(b)-Dl(d) and are given in figure D2(a) 
for the path axis components of the apparent dispersions (6Rp,6V p ,6a p ) ; for the true 
dispersions, (6R p , 6V p , 6a p ) ; and for the control activity, including the net activity, 
Cp = Aa p + b p (Coriolis effects are neglected and results apply to the linear (unsat- 
urated) control domain). The steady state step responses are included in the table 
and these responses show that constant estimation errors for any state, R,V,a, are 
converted by the control to a position offset, SR, with no steady state velocity 
error or net control offset. In effect, the control can detect and null out constant 
velocity and acceleration estimation errors at the cost of a position hang off and 
it leaves information on these biases in the steady state apparent dispersions and 
control offsets. 

Frequency response magnitudes provide further information on excitation of the 
system by the estimation errors. Analytical expressions are given in figure D2(b) 
for the true dispersions and net control activity; as seen, all output variables 
depend on all input estimation errors (6R depends on R,V,a, and similarly for 6V,c) 
and the asymptotic behavior (u) -*■ <») indicates nonzero control activity in response to 
R,V over the entire frequency spectrum. Numerical amplitude ratio results are shown 
in figure D2(c) for several values of the control bandwidth. These results indicate 
the following trends. First, steady state position tracking error responds 


124 



APPARENT 

DISPERSIONS 


TRUE 

DISPERSIONS 


CONTROLS 


6R 

A 

6V 

•A 

6a 

6R 

5V 

6a 

Aa 

b 




A, 

q 

Kv 

i 

q 

q(rs+ 1) 

a 

M 

* L* 
1 

s 

s — 

q 

q(rs + 1) 

0 

0 

rs 

rs + 1 

A 

A 

1 

q 

q 

q(rs + 1) 

S A 

S A 

s 

q 

q 

q(rs + 1 ) 

s 2 A 

s 2^ 

$2 

q 

q 

q(rs + 1) 

. s 2A 

$2 A 

1'4)^ 

q 

q 

0 

0 

rs + 1 

2 A 

s 2 A 

s 2 

-5 

q 

-s — 
q 

q(7S + 1 ) 


TRANSFER FUNCTIONS 


NOTES: 

q = s 2 + sKy + 
6R = RC - R 
6V = VC- V 
6a = act - a 
c = Aa + b 


0 

K v 

K r 

1 

K r 

0 

-1 

0 

0 

0 

0 

1 

K v 

1 

K r 

K r 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

-1 

0 

0 

0 

. 


_ 


R 

V 

a 


STEADY STATE STEP RESPONSE 


(a) Linear domain transfer functions. 
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(b) Frequency-response amplitude ratios. 


Figure D2.- Trajectory-t racking errors and control activity resulting from 

estimation errors. 
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(c) Frequency-response amplitude ratios: p , t = 0.707, 5 sec. 

Figure D2.- Concluded. 

principally to low-frequency estimation errors (below w n ) , with maximum sensitivity 
at o) ft = 0. Bias errors in any state estimate are mapped by the control into posi- 
tion hangoff, and higher-frequency errors are attenuated so that position dispersions 
will be "smoother" than the estimation errors. Second, steady state velocity track- 
ing errors are unaffected by estimation error biases and have maximum sensitivity to 
errors at the control bandwidth (w = w n ) ; sensitivity to V is 1.0 at w n and 
otherwise velocity dispersions reproduce principally the estimation error content at 
frequencies near o) n . And third, the net control activity is insensitive to low- 
frequency errors; it responds principally to frequencies at and above the control 
bandwidth and reflects the higher frequency content of R and V. 

Thus, it is the frequency of the estimation error relative to the control band- 
width that determines its effect on the system outputs. In this regard, we note that 
flight safety considerations depend principally on position tracking accuracy and, 
therefore, on the low-frequency estimation errors; on the other hand, velocity 
variations, control activity, and ride quality depend principally on the estimation 
error spectrum at and above o) n . 
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SIMULATION RESULTS: TRACKING ERRORS 


Ensemble rms and sample case trajectory tracking errors excited by the navigation 
errors are given in figures D3 and D4, along with a comparison with the corresponding 
state estimation errors. These results contain various effects not studied in the 
response analysis described above, including control system transient responses, con- 
trol saturation, neglected Coriolis effects, and the collective effects of mutually 
correlated state estimation errors. 

The position dispersions show several features. First, rms(6R) exceeds rms(R) 
significantly and persistently in response to large changes in rms(R), such as the 
lateral axis ramp reduction during leg 1 and the step changes on all axes at entry to 
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Figure D3.- System response to estimation errors: ensemble rms, 
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Figure D4.- System response to estimation errors: sample case. 


MODILS or elevation coverage. During these periods, the control^system maneuvers 
within the rate limit, to null the apparent dispersions, <$R, so that 6R 

necessarily lags R. Excursions in excess of rms(R) are also excited on the longi- 
tudinal and lateral axes during turns because of (1) low-frequency changes in R 
imposed by rotation of the path axes and (2) cross-coupling of the^path axes disper- 
sions by uncompensated Coriolis effects. Second, rms(6R) and rms(R) are the same 
size where the estimation error is predominantly bias (examples are seen in both rms 
and sample case results for the normal and longitudinal axes before calibration of 
the baroaltimeter and VORTAC range biases, respectively). This result also occurs on 
the glide slope where R variations are low in amplitude and 6R reproduces the low- 
frequency content of R with some lag and little amplification (e.g. , the sample case 
normal axis errors). Third, the sample case history shows that SRp(t) is signifi- 
cantly smoother than Rp(t), as expected. 
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The velocity dispersion satisfies the identity 


6V = -V + 6V 
P P P 

and reproduces^ the estimation error frequency content near the control bandwidth 
(principally V in the frequency range 0.5 w n to 3 w n ) combined with the transient 
maneuvering velocity needed to null the apparent position dispersion. The control law 
limits the maneuvering dispersion to SVl^. 

The rms velocity excursions show excursions well in excess of rms(V) during 
large changes in rms(R) (e.g., on the longitudinal axis when passing near the 
VORTAC and on all axes at entry to MODILS or elevation coverage). These excursions 
reflect the maneuvering needed to null the apparent position dispersion. The rms 
lateral velocity dispersion is initially smaller than rms(V) as a result of the con- 
trol law’s ability to null the lateral velocity estimation^bias owing to bearing bias. 
More generally, rms(6V) shows persistent excess over rms(V) by a factor of 1.5 to 
2.0, because the velocity estimate is generally unbiased, except as noted above, and 
because the control is continually excited by estimation errors. The sample case 
illustrates the underlying behavior in detail; the normal axis velocity dispersion is 
seen to reproduce V with some amplification of most peaks during the extended 
periods of stationary errors, and a low-frequency excursion is superposed on this 
behavior at entry to MODILS coverage. The longitudinal and lateral axis results show 
additional examples of the behavior described above. 

Ensemble extremes of position and velocity dispersions were reviewed and found 
generally to be twice the rms dispersions given in figure D3. These extreme disper- 
sions are due to the estimation errors of the present system and are experienced on 
one or more approaches out of 10 if they depend principally on sample case constants, 
or several times during an approach if they depend on random measurement noise. 


SIMULATION RESULTS: CONTROL ACTIVITY 


The net control activity excited by estimation errors is illustrated by the rms 
and sample-case results shown in figure D5 . In conventional aircraft, these histories 
map into corresponding engine, roll angle, and pitch angle activity. Control activity 
depends on sensor accuracy, on the control bandwidth, and on details of the estimation 
errors. The rms longitudinal control for the present 4D guidance law is roughly uni- 
form at an average level of 0.04 g throughout the approach; the normal axis, with 
better position and acceleration estimation accuracy, shows less activity, at about 
0.02 g (this difference would be greater if the control bandwidth were identical for 
the two axes). On the other hand, longitudinal control activity is about the same 
during both MODILS and VORTAC use, indicating an insensitivity to position accuracy in 
the test context. To clarify the dependence on sensor accuracy, results for a high- 
accuracy IMU are included in figure D5. When results for the two IMU’s during MODILS 
use are compared, it is apparent that control activity depends strongly on IMU accuracy 
when position accuracy is good. Similarly, a comparison of results during VORTAC use 
with those during MODILS use for the inertial grade IMU indicates a strong dependence 
on position estimation accuracy when IMU accuracy is good. Thus, for the present sys- 
tem, position sensor accuracy dominates longitudinal and lateral control requirements 
during VORTAC use, and the IMU accuracy dominates it during MODILS use. Finally, it 
is apparent that both the VORTAC and IMU accuracy must be improved to obtain signifi- 
cant reduction in control activity because of estimation errors throughout the termi- 
nal area. 
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Figure D5.- System response to estimation errors: net control activity 

({a) n .} = (0.125, 0.125, 0.3} ). 

The regulator control authority is typically 0.1 g; the rms control activity in 
figure D5 indicates that 20% to 40% of this authority, depending on control axis, is 
used up by the estimation errors for the present system but that this could be 
reduced to negligible levels compared to the authority with improved sensors. For 
aircraft, turbulence maps into acceleration disturbances oriented principally along 
the normal axis so that a sufficient reserve of control margin for external distur- 
bances is more important for this axis. For the remaining axes, it is likely that 
estimation errors excite a larger average disturbance, Aa(R,V,a), than does the tur- 
bulence field. Control activity that results from estimation errors also depends on 
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the control bandwidth o) n and increases 
with o) n , as does the frequency of satura- 
tion, while reserve margin and dispersion 
peaks decrease. Some results for the 
effect of bandwidth on rms control activ- 
ity (averaged over the entire approach) 
are given in figure D6 for the present 
system; it illustrates this trend and also 
confirms the distinctly lower control 
activity for the normal axis at a given 
w n because of its distinctly better esti- 
mation error environment. 

The sample case history (fig. D5) 
shows typical details: the dominant fre- 

quencies are seen to be higher than those 
of the dispersions 6Rp,6V p , as expected, 
and saturation of the control authority 
occurs infrequently on the longitudinal 
and lateral axes and not at all on the 
normal axis. 


Ui 


>- 
C/5 h- 



TRANSIENT RESPONSE NATURAL 
FREQUENCY - rps 

Figure D6.- Average control activity. 


EFFECT OF ATTITUDE ESTIMATION ERRORS ON SYSTEM PERFORMANCE 


A detailed study of the effects of inner loop sensor and estimation errors on 
system performance is beyond the scope of our simulation, but some general comments 
can be made. 

Inner loop estimation errors result in acceleration control errors, f, which 
induce dispersions and compensation activity by the translational control in accor- 
dance with the model of figure Dl(d). In particular, the attitude measurements pro- 
vided by the gyroscopes of the present system have low-frequency errors that can be 
assumed to result in like attitude control errors and corresponding errors in con- 
trolling aerodynamic force. These latter errors induce dispersions and control 
activity, depending on frequency, but can be compensated by the plant model error 
compensator l>, up to the limit of its authority (0.1 to 0.2 g) . An important effect 
of this compensation activity is the corresponding loss of margin to compensate for 
aircraft identification errors in the plant inverse portion of the control laws, 

G(z ,a) . 

Roll angle determines lift force orientation about the longitudinal axis and is 
the primary control over lateral accelerations; roll errors of 0.5° result in accel- 
eration errors of 0.01 g — an appropriate accuracy level for the roll control. In 
the present tests, the roll angle measurement shows poorer accuracy, with excursions 
to ±1.5° during turns. 

Angle of attack, which is controlled through pitch attitude, determines the lift 
force magnitude and is the primary control over normal acceleration in conventional 
aircraft. Lift sensitivity to angle-of-attack control errors is C^/Cl g/rad. 

This ratio typically increases with airspeed and was found to be in the range 0.04 to 
0.2 g/deg over the speed range of the test approach path for an example STOL aircraft. 
For this example, angle-of-attack estimation accuracies of 0.25° to 0.05° are needed, 
depending on airspeed, to obtain lift control accuracy at 0.01 g. Angle-of-attack 
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error is a superposition of several errors (eq. (5.12)), with pitch measurement errors 
dominating the low-frequency content. In the present test, pitch measurement error 
showed excursions over ±1.5° during maneuvering; the excursions were independent of 
speed, and are, therefore, an important outer loop disturbance, particularly in 
higher speed maneuvering where the error is an order of magnitude larger than the 
desired accuracy. 

Last, sideslip angle determines side force and is controlled through heading 
angle. However, sideslip is unobservable to the present estimation system and cannot 
be controlled. In this case, the automatic control relies on the aircraft’s weather- 
vane stability to null side force in steady state, and the control is insensitive to 
the large error in the directional gyroscope. 
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